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ABSTRACT 


Digital  control  offers  exciting  possibilities  for  ac¬ 
curacy  and  flexibility  that  cannot  be  achieved  by  analog 
means.  The  USQ-2Q  General  Purpose  Digital  Computer  is  cn 
board  several  ships  today  for  the  Navy  Tactical  Data  System, 
and  an  effort  is  made  to  show  how  such  a  computer  might  also 
be  employed  directly  in  closed-loop  shipboard  fixe  control 
systems.  Some  of  the  theory  of  sampled  data  systems  and 
digital  compensation  is  discussed.  A  digital  computer  pro¬ 
gram  which  simulates  the  action  of  a  sampled  data  system 
with  a  digital  controller  is  employed  to  experimentally  de¬ 
termine  the  behaviour  of  various  physical  plants  subject  to 
digital  control.  An  evaluation  is  made  of  the  noise  rejec¬ 
tion  capabilities,  steady  state  error,  and  transient  control 
of  a  digital  process.  A  full  description  and  instructions 
for  its  use  is  given  for  the  simulator  program,  which  has 
proved  to  be  very  helpful  in  making  theoretical  studies  of 
control  system  behaviour. 

The  writer  wishes  to  express  his  appreciation  for  the 
assistance  and  encouragement  given  him  by  Professor  Mitchell 
L.  Cotton  of  the  U„  S„  Naval  Postgraduate  School  and  John  B. 
Slaughter  of  the  U„  Sc  Navy  Electronics  Laboratory  in  this 
investigation. 
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Introduction 


A  control  system  which  contains  a  digital  computer 
may  be  classed  as  a  sampled-data  system*  Such  systems  are 
becoming  increasingly  more  common,  both  in  military  opera¬ 
tions  and  in  the  commercial  field.  By  the  use  of  sampling 
in  control  systems  it  is  possible  to  build  simple,  sensi¬ 
tive,  and  efficient  power  control  devices.  Sampling  per¬ 
mits  the  control  of  tremendous  power  by  sensitive  control 
elements  without  excessive  amplification,  and  also  mini¬ 
mizes  the  loading  effect  upon  sensitive  instruments. 
Utilizing  sampled  data  and  digital  components  in  a  control 
system  allows  time  sharing  of  important  parts  of  the  system. 
This  is  a  very  significant  advantage  of  digital  control 
systems. '/j'/ 

Since  sampled-data  can  readily  be  coded,  the  data 
signals  in  digital  control  systems  are  received  and  trans¬ 
mitted  in  pulse-code  form,  and  they  will  provide  almost 
error-free  channels  for  transmission  through  noisy  media. 
Consequently,  the  only  noise  in  the  transmission  of  pulse- 
coded  signals  is  the  quantization  error.  Likewise,  a  digi¬ 
tal  computer  in  a  control  system  allows  data  processing  of 
control  information  so  that  the  flexibility  and  versatility 
of  the  digital  computer  can  be  utilized  to  improve  the  per¬ 
formance  of  the  control  system.  Furthermore,  digital  con¬ 
trol  techniques  make  feasible  the  system  compensation  by 
non-linear  programming  and  adaptive,  or  self -optimizing, 
control.  This  is  another  significant  advantage  achievable 
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with  digital  control 


/2/ 


The  principal  features  of  a  typical  system  containing 


a  digital  computer  are  shown  in  block  diagram  form  in  Fig- 
ure  1 . 
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Figure  1.  Typical  block  diagram  of  a  hybrid  control  system. 

A  characteristic  feature  of  this  system  is  the  fact  that 
data  will  appear  at  one  or  more  places  within  it  as  a  se¬ 
quence  of  numbers  or  fixed  signal  levels  that  do  not  change 
appreciably  except  at  fixed  instants  in  time.  The  process 
that  converts  continuous  data  into  such  a  sequence  is  called 
the  sampling  process,  and  the  mechanism  that  implements  the 
process  is  called  an  analog- to-digital  converter  (ADC).  The 
process  that  converts  the  output  number  sequence  of  the  digi¬ 
tal  computer  to  a  more  or  less  continuous  analog  quantity 
is  called  the  desampling  process,  and  the  mechanism  accomp¬ 
lishing  this  process  is  called  the  digital-to-analog  con¬ 
verter  (dac  )  y2^ 

The  quantities  T^,  T2,  and  T^  are  the  periods  of  the 
samplers  shown.  In  almost  every  case  T^  is  the  same  as  T2, 
but  the  relation  between  T^  and  T^  may  be  chosen  by  the  de¬ 
signer.  The  simplest  relation,  of  course,  is  when  T?  equals 
T^.  The  system  is  then  called  single-rate  since  all  discrete 
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number  sequences  appearing  within  the  system  have  the  same 


sampling  period.  A  multi-rate  system  occurs  when  T?  and 
are  not  the  same.  Analysis  is  simplified  when  and 
are  related  by  an  integer  and  we  have  11  synchronous" 
sampling.  It  is  also  simpler  to  analyze  a  multi-rate  sys¬ 
tem  wherein  T3  is  less  than  T 2  (ie.,  T3  -  T^/n  n  =  2,3,...) 
which  is  fortunate,  since  most  of  the  advantages  of  a  multi¬ 
rate  system  are  derived  when  the  output  of  the  digital  con¬ 
troller  occurs  at  a  faster  rate  than  the  input. 

The  experimental  results  described  in  this  thesis  are 
obtained  from  a  digital  computer  program  that  simulates  a 
hybrid  control  system  from  beginning  to  end.  See  Appendix  I 
for  the  description  of  this  program.  In  an  actual  control 
system  the  only  operations  that  would  perforce  be  external 
to  the  digital  computer  would  be  the  plant  itself  and  the 
ADC.  In  the  investigation  conducted  by  John  Slaughter  at 
the  U.  S.  Navy  Electronics  Laboratory,  San  Diego,  California, 
a  Packard-Bell  Model  M2  analog- to-digital  converter  is  em¬ 
ployed,  and  a  100  micro-second  conversion  and  settling  time 
is  allowed.  For  a  discussion  of  a  proposed  analog- to-digital 
system  for  use  with  the  CDC  1604  digital  computer  see  Ref¬ 
erence  6. 

The  generation  of  the  error  signal  by  comparison  be¬ 
tween  input  and  output  signals  can  be  accomplished  either 
externally  to  the  digital  computer  or  by  the  digital  com¬ 
puter  itself.  The  latter  case  would  require  a  two  channel 
ADC  and  a  multiplexing  system  for  sampling  respectively  the 
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input  and  the  output  signals. 

Data  reconstruction  by  the  DAC  can  also  be  performed 
by  an  external  device  and  can  be  of  various  degrees  of  com- 
plexity.  However,  the  fact  that  the  plant  in  control  systems 
is  usually  low-pass  in  frequency  response  makes  the  use  of 
overly  complex  data-reconstruction  systems  unnecessary. 
Generally  speaking,  the  zero-order  data  hold  is  found  to  be 
adequate  for  most  systems  found  in  practice.  This  being  so, 
the  data  reconstruction  process  can  easily  be  implemented 
by  the  digital  computer  since  an  output  register  of  the  com¬ 
puter  can  be  utilized  as  a  zero-order  hold. 
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2.  The  effect  of  sampling. 

Sampling  is  the  act  of  examining  a  continuous  function 
at  specified  increments  of  an  independent  variable,  usually 
time.  Sampling  processes  can  and  do  assume  a  wide  variety 
of  forms  and  can  be  separated  roughly  into  three  categories 

1.  Linear,  fixed  pattern  sampling 

2.  Non-linear,  signal  dependent  sampling 

3.  Random  sampling 

Included  in  the  first  category  is  the  historically 
conventional  concept  of  "instantaneous"  sampling  of  fixed 
period  T.  By  instantaneous  it  is  meant  that  the  duration, 
or  pulse  width,  of  the  sample  of  the  signal  is  very  short 
compared  to  the  shortest  time  constant  of  the  system*  Un¬ 
der  this  assumption  (which  will  be  employed  in  this  thesis) 
the  sampling  process  can  be  regarded  as  generating  a  se¬ 
quence  of  impulses.  This  viewpoint  leads  to  the  develop¬ 
ment  of  the  Z-transform  calculus  and  also  to  the  use  of 
difference  equations  as  a  mathematical  tool.  Extensive 
tables  of  Z-transforms  may  be  found  in  the  literature,  and 
it  is  relatively  easy  to  take  the  Z-transform  of  simple, 
low-order  transfer  functions.  However,  for  a  fifth  order 
plant  or  higher,  finding  the  Z-transform  is  quite  laborious 
Appendix  2  describes  a  digital  computer  program  written  for 
transfer  functions  of  Type  0  and  Type  1  servo-mechanisms 
up  to  the  tenth  order. 

The  Sampling  Theorem  by  Shannon  states  that  a  band- 
limited  signal  of  highest  frequency  Fq  can  be  completely 
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specified  by  sampling  at  a  rate  of  2Fq.  Recovery  of  the 
signal  at  this  minimum  theoretical  sampling  rate  requires 
an  ideal  filter.  In  practice,  signals  are  not  ordinarily 
band-limited  and  filters  do  not  have  ideal  attenuation 
characteristics.  Even  at  high  sampling  rates,  sampled- 
data  systems  usually  exhibit  distortion  in  the  form  of 
"ripple".  One  means  of  combating  this  ripple  is  to  employ 
a  multi-rate  digital  controller  whose  output  occurs  at  a 
faster  rate  than  the  basic  error-sampling  rate  of  the  over¬ 
all  system.  By  applying  control  signals  to  the  plant  more 
frequently  much  of  the  ripple  can  be  reduced  and  the  plant 
can  be  made  to  respond  more  quickly  to  an  input  signal. 
Design  judgment  must  be  used,  however,  since  the  output 
signals  from  a  multi-rate  controller  are  frequently  of 
greater  magnitude  than  those  from  a  single-rate  device, 
and  the  possibility  of  plant  saturation  as  well  as  the 
added  complexity  of  the  controller  must  be  balanced  against 
the  advantages  to  be  gained  from  a  multi-rate  system. 

Of  crucial  importance  in  a  sampled-data  system  is  the 

effect  of  sampling  on  the  stability  of  the  system.  In  an 

otherwise  stable  system,  sampling  reduces  the  degree  of 

/2/ 

stability  and  can  even  lead  to  instability.  '  There  are 
various  criteria  for  continuous  systems  that  can  be  applied 
to  sampled-data  systems.  Among  these  are  the  modified 
Routh-Hurwitz  criterion,  the  Schur-Cohn  criterion,  and  the 
Nyquist  criterion.  As  an  example  of  the  determination  of 
stability  using  the  modified  Routh-Hurwitz  criterion  we 

will  use  the  third  order  Type  1  servo-mechanism  shown  in 
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Figure  2.  The  plant  is  closed-loop  stable  with  the  gain 
indicated.  We  desire  to  convert  the  continuous  system  to 
a  sarnpled-data  system/  so  we  first  insert  sampling  and  data 
reconstruction  as  shown.  The  matter  of  digital  compensation 
will  be  held  in  abeyance  while  we  investigate  the  effect  of 
sampling  only. 
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Figure  2.  A  third  order  sampled-data  system. 


Suppose  that  a  sampling  rate  of  100  samples  per  sec¬ 
ond,  or  T  .01  seconds,  has  been  chosen  a  priori.  Then 
the  modified  Routh-Hurwitz  procedure  gives  a  "yes  or  no" 
result  in  regard  to  stability.  For  this  value  of  T  the 
forward  path  pulse  transfer  function  (or  Z-transform)  is: 

P(-a)  .0 IO<Z,OOBZ  +-  .033672.  2  +  ^00  6^^4 

Z3-  3. 31S38&Z*'  +  /.  £6  2/27  2  -O. 

For  the  unity  feedback  configuration  shown  the  closed  loop 

transfer  function  has  a  characteristic  equation  given  by: 

P(z)  +  Q( z)  =  A(z)  =  A3  z^  +  A2  z^  +  Al  z  +  A0 

where  A3  =  1.000000 

A2  =  -2.264786 
Al  =  1.695799 

A0  =  -0.380157 

The  closed  loop  sampled-data  system  is  stable  if  the  roots 
of  A(z)  lie  within  the  unit  circle  in  the  Z-plane.  By 
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employing  the  Bilinear  Transformation,  Z  =  1  +  w,  the  in- 

1  -  w 

terior  of  the  Z-plane  is  mapped  into  the  left  half  of  the 

w-plane.  Once  the  modified  characteristics  equation  is 

3  2 

obtained,  namely:  A(w)  =  B3  w  +  B2  w  +  B1  w  +  BO  then 

the  Routh  array  may  be  formed  and  the  criterion  for  a  third 

order  system  may  be  applied.  This  criterion  is: 

Bl*B2  -  B0*B3  is  a  positive  quantity  if  all  the  roots 
of  A(w)  are  in  the  left  half  plane 

In  this  example,  B3  =  5,340742 

B2  =  2.428516 
B1  =  0.179886 
BO  =  0.050856 

so  B1*B2  -  30*B3  =  +  0.165247  and  the  system  is  stable  for 
T  equals  .01  seconds. 

If  it  is  desired  to  decrease  the  sampling  rate  as  much 
as  possible  in  order  to  time- share  the  computer  or  to  allow 
for  increased  computation  time  between  samples,  then  it  be¬ 
comes  necessary  to  know  the  limit  of  stability  as  the  samp¬ 
ling  period  T  is  increased.  Unfortunately,  the  character¬ 
istic  equation  is  transcendental  in  T,  and  we  cannot  solve 
directly  for  T  at  the  stability  limit.  Iterative  procedures 
most  readily  delegated  to  a  digital  computer  must  be  em¬ 
ployed.  Appendix  3  shows  the  development  of  the  character¬ 
istic  equation  and  the  Fortran  program  for  the  CDC  1604 
digital  computer  which  calculates  the  Routh-Hurwitz  criter¬ 
ion  for  this  third  order  system  for  one  hundred  values  of  T 
ranging  from  ,001  seconds  to  0.1  seconds.  From  the  results 
shown  the  stability  limit  on  T  is  evidently  near  .025  sec¬ 
onds,  since  this  is  the  last  value  for  which  the  R-H  criter¬ 
ion  is  a  positive  quantity. 
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Figure  3  shows  the  response  of  this  un-compen sated 
third  order  sampled  data  system  to  a  unit  step  input  for 
three  values  of  T.  For  T  of  .001  seconds  the  response  is 
indistinguishable  to  the  eye  from  the  results  obtained  in 
an  analog  computer  simulation  of  the  problem.  This  would 
be  a  good  value  to  choose  for  the  time  increment  in  a 
solution  of  the  plant  differential  equation  by  numerical 
methods  (see  Appendix  1).  It  can  be  seen  that  the  sta¬ 
bility  of  the  system  is  degraded  somewhat  when  T  is  in¬ 
creased  to  .01  seconds,  but  the  system  is  still  stable. 
This  is  the  sampling  period  to  be  employed  later  when  we 
introduce  digital  compensation  to  this  third  order  system. 
A  sampling  period  of  .05  seconds  results  in  severe  insta¬ 
bility,  as  we  could  expect  from  the  Routh-Hurwitz  deter¬ 
mination.  Digital  compensation  is  a  must  in  this  case, 
since  stability  is  a  prime  requisite.  Whether  or  not  we 
compensate  a  stable  system  depends  on  how  satisfactory  the 
response  is  without  compensation. 


9 


10 


3.  The  digital  process. 

As  previously  mentioned,  a  sampled-data  system  is 
characterized  by  the  fact  that  in  certain  areas  of  the  sys¬ 
tem  signals  appear  as  discrete  quantities  or  fixed  numbers 
in  a  sequence  at  regular  intervals  of  time.  The  most  ob¬ 
vious  point  at  which  this  occurs  is  at  the  digital  computer 
or  controller.  The  input  to  the  computer  is  a  sequence  of 
numbers  and  the  computer  output  is  also  a  sequence  of  num¬ 
bers.  We  can  thus  regard  the  computer  as  a  "black  box"  as 
shown  in  Figure  4. 


ElOTt] 

Ducita  l 

E2[pT2] 

Computer 

p  =  1,  t,  . ,  M 

Figure  4.  Generalized  concept  of  a  digital  computer. 

The  input  sequence,  El  0TJ  •  consists  of  N  discrete 
values  separated  in  time  by  T^  seconds.  The  output  se¬ 
quence  may  contain  a  different  number,  M,  of  fixed  quan¬ 
tities  and  the  period,  T2,  of  the  output  sequence  need  not 
be  the  same  as  the  input  period. 

For  simplicity  let  us  first  take  T,  equal  to  T^. 

Then  at  a  certain  time,  pT,  the  general  form  of  linear 
computation  that  the  digital  computer  can  perform  may  be 
expressed  ass 

e2 i>T]  +  bi-e2[(f>-0T]  +•••+-  bM_( •  el  [Cp-w+i)t] 

-  clQ  •  el  [pf]  t  ai‘  el  [(P-OT]  £(M-ei[(p-N  +  «)T] 
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This  can  be  written: 


A 1-1 


el  OtJ  =  2_  etf-K 

M-  o 


M-i 


Z  bJ  e2p-j 

J=< 


or,  by  Z-transform  calculus: 


^Z(s) 

£1(?) 


Z^kZ 

1  +  Zbj  *  J 


which  can  be  considered  the  pulse  transfer  function  of  a 
digital  computer.  If  all  the  b^  are  zero  then  we  have  what 
is  called  a  finite  memory  process  since  the  output  at  pre¬ 
sent  time,  pT,  is  the  weighted  sum  of  the  present  input 
plus  earlier  inputs  up  to  some  finite  limit,  N.  If  not  all 
of  the  b  coefficients  are  zero  then  a  weighted  contribution 
of  earlier  outputs  is  also  involved  in  forming  the  output 
at  time  pT.  Since  each  earlier  output  was  a  combination  of 
both  input  and  output  sequences,  ad  infinitum,  we  now  have 
what  is  called  an  infinite  memory  process.  Another  way  of 
showing  this  is  to  perform  long  division  on  the  transfer 

function  to  obtain  new  weighting  coefficients,  c^.,  for  the 

_ 

input  sequence  alone.  Only  in  the  case  where  a^Z  con¬ 
tains  1  -  bjZ”-1  as  a  factor  does  long  division  fail  to 
result  in  an  infinite  series,  showing  the  need  for  an  in¬ 
finite  memory  of  previous  inputs.  When  such  a  transfer 
function  describes  the  over-all  closed  loop  process,  we 
have  what  is  known  as  a  "non-finite  settling  time"  process 
and  the  systematic  error  (or  difference  between  input  and 
output)  can  never  theoretically  reach  zero. 

In  the  system  configuration  employed  in  this  thesis 
we  will  use  the  concept  of  the  digital  computer  transfer 
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function  in  two  ways.  First  is  the  over-all  pulse  transfer 
function  of  the  closed-loop  process,,  defined  ass 

/?(&)  ~  //JPUT  SSOUSTAJC-e 
Our  attention  will  be  devoted  mostly  to  finite  memory  K(z) 

consisting  of  numerator  polynomial  only.  Such  a  process 
has  a  finite  settling  time  and  can  "forget"  the  effect  of 
initial  errors  or  acquisition  transients  in  the  input  sig¬ 
nal. 


Second  is  the  pulse  transfer  function  of  the  digital 
controller  in  the  forward  path.  Storage  will  be  allocated 
for  a  finite  number  of  values  in  both  the  input  error  se¬ 
quence  and  the  output  error  sequence,  and  the  controller 

-1 

transfer  function  will  be  a  ratio  of  polynomials  in  Z 

defined  as:  _ ,  x  _  Eout(z)  _  processed  error  sequence 
'z'  Ein( z)  observed  error  sequence 

In  case  the  input  sequence  to  the  controller  has  a 
period,  T^,  different  from  the  period,  T2,  of  the  output 
sequence  we  can  still  write  a  pulse  transfer  function  that 
is  amenable  to  interpretation  using  difference  equations 
providing  and  T2  are  related  by  an  integer.  The  pulse 
transfer  function  is  then  expressed  in  the  "Z-domain" 
created  by  the  faster  sampler,  which  in  this  thesis  will 
be  at  the  output  from  the  controller.  In  the  next  section 
is  provided  the  mathematical  basis  for  these  assertions. 
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4.  Response  criteria. 

In  view  of  the  flexibility  possible  with  active  digi¬ 
tal  compensation  there  is  a  very  large  number  of  possible 
over-all  response  functions  which  can  be  implemented.  As 
a  starting  point,  the  simplest  or  "minimal "  prototype  re¬ 
sponse  functions  are  convenient.  Minimal  prototype  systems 
are  approached  with  the  point  of  view  that  they  must  be 
able  to  respond  satisfactorily  to  some  convenient  test  in¬ 
put  such  as  a  step,  ramp,  or  constant  acceleration,  or  all 
three.  The  requirements  which  are  set  for  minimal  prototype 
response  functions  are: 

a.  The  over-all  response  and  the  response  of  all  ele¬ 
ments  of  the  system  must  be  physically  realizable. 

b.  The  steady-state  response  to  the  test  input  must 
have  zero  systematic  error. 

c.  The  transient  response  should  be  as  fast  as  pos¬ 
sible  and  the  settling  time  should  be  equal  to  a 

/4/ 

finite  number  of  sampling  intervals. 

The  following  table  shows  the  test  inputs  and  the  required 
minimal  prototype  response  functions,  K(z): 


Input 

K(z) 

Step 

-1 

z 

^  -1 

-2 

Ramp 

2z 

-  z 

Acceleration 

3z 

,  -2  ,  -3 

-  3z  +  z 

These  minimal  responses  are  possible  only  for  those  plants 
whose  pulse  transfer  functions,  G(z),  contain  no  zeros  or 
poles  outside  the  unit  circle  in  the  2-plane. 
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The  essence  of  digital  compensation  is  that  the  con* 
stant  coefficients  of  the  pulse  transfer  function,  D(z), 
for  the  controller  are  adjusted  so  that  the  desired  over* 
all  pulse  transfer  function  of  the  process,  K(z),  is  real¬ 
ized.  Depending  on  the  structure  of  the  sampled-data  feed¬ 
back  system,  the  controller  transfer  function  is  derived 
from  system  functions  in  various  ways.  The  configuration 
employed  in  this  thesis  is  the  typical  one  shown  earlier  in 
Figure  1,  viz.,  in-line  digital  compensation  of  a  unity 
feed-back,  error-sampled  control  system.  The  ADC  can  be 
thought  of  as  a  single  sampling  switch  operating  with  per¬ 
iod  T.  Data  reconstruction  will  be  accomplished  by  a  zero- 
order  hold,  and  the  pulse  transfer  function  of  this  DAC  in 
combination  with  the  continuous  plant  will  be  denoted  by 
G(z).  Since  the  actual  input  and  output  are  continuous 
signals,  fictitious  samplers  are  shown  producing  R(z)  and 
C(z)  so  that  we  may  speak  of  the  over-all  pulse  transfer 
function,  K(z).  For  convenience.  Figure  1  is  re-drawn 
embodying  these  ideas  as  Figure  5  below. 


o— 


Digital 


PlAi^t  (iNCLuoes 


V  E1(Z> 

P(z) 

\  E2  (a) 

?-  T 

r/n 

Gr 

t 

_ J 

K(^) 


cC-fc) 

-o 


,  CCz) 

- o 


Figure  5.  Sampled-data  system  with  in-line  compensation. 


Taking  the  single-rate  case  first  (where  n-1,  thus  T/n-T), 
the  over-all  pulse  transfer  function  is  seen  to  be 
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K(z) 


C(  z)  D(z)«G(z) 

R(zT  1  +  d(z)«g(z) 


which  can  be  solved  to  give 


D(z) 


1  K(z) 

gTTT  *  i-,”k(z7 


Depending  on  the  poles  and  zeros  of  the  plant  which 

it  is  desired  to  compensate  digitally*  complete  freedom  of 

choice  of  K(z)  may  not  exist.  According  to  Franklin  and 

Ragazzini  (4)  the  restrictions  on  K(z)  ares 

It  is  necessary  that  the  specified  over-all  pulse 
transfer  function,  K(z),  contain  as  its  zeros  all 
those  zeros  of  the  plant  pulse  transfer  function, 
G(z)  which  lie  on  or  outside  the  unit  circle  in 
the  z-plane,  and  that  l-K(z)  contain  as  its  zeros 
all  those  poles  of  the  plant  transfer  function  which 
lie  on  or  outside  the  unit  circle  in  the  z=plane. 

Systems  which  are  designed  for  minimum  and  finite 
settling  time  often  do  not  give  good  performance  when  sub¬ 
jected  to  an  input  other  than  that  for  which  they  are  de~ 
signed.  In  this  sense,  minimal  systems  may  be  regarded 
as  highly  "tuned".  In  addition,  the  severe  shocks  which 
result  in  the  plant  cause  substantial  ripple  in  the  con¬ 
tinuous  output  even  though  the  output  is  error-free  at  the 

sampling  instants.  These  effects  were  noted  in  the  lit- 
/7  / 

erature  and  led  to  the  introduction  of  a  term  in  the 
over-all  response  function  known  as  the  staleness  factor. 
This  factor  leads  to  a  softening  of  the  response,  with  the 
result  that  the  system  can  be  expected  to  respond  slightly 
slower  but  adequately  to  a  number  of  test  inputs.  The 
choice  of  staleness  factor  can  be  arrived  at  by  optimizing 

/g  / 

procedures'  '  or  by  observing  the  response  to  some  most 
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likely  form  of  test  input.  The  staleness  factor  is  intro¬ 
duced  as  the  constant  "c"  in  the  following  expression  for 
the  over-all  prototype  transfer  function: 

tw_\  _  a»K(  z)rainimal 
Z'  ~  —IN 

(l-cz-V 

where  in  practice,  N  =  1,  c  is  a  positive  constant  less 
than  unity,  and  "a"  is  determined  by  applying  the  Final 
Value  Theorem. 

In  order  to  compensate  the  system  so  that  the  output 
is  ripple-free  after  a  finite  settling  time  the  following 

/4/ 

design  rules  must  be  applied;' 

a.  All  the  rules  for  minimal  prototype  response 
systems  apply. 

b.  The  open  loop  transfer  function  must  be  capable 
of  generating  an  output  that  is  the  same  as  the 
input. 

c.  The  over-all  pulse  transfer  function,  K(z),  must 
contain  as  its  zeros  all  the  zeros  of  the  plant 
pulse  transfer  function,  G(z),  and  not  just  the 
zeros  of  G(z)  which  lie  outside  the  unit  circle 
in  the  z-plane. 

Let  us  turn  now  to  the  case  in  which  the  output  of 
the  digital  controller  occurs  more  frequently  than  the  in¬ 
put.  The  development  that  follows  is  due  to  Kranc  (Refer¬ 
ence  9 ) . 

If  we  apply  a  discrete  sequence,  E(z),  to  a  continuous 
plant,  H(s),  then  it  follows  that  the  continuous  output  is; 

C(s)  =  H(s)-S(z)  (1) 

For  n  equal  unity  we  may  take  the  pulse  transform  of  both 
sides  of  equation  (1)  for  a  sample  period  of  T.  Since  E(z) 
is  already  a  function  of  Z  it  is  invariant  under  the  trans- 
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formation,  so  the  result  iss  C(z)  =  H(z)»E(z) 

When  the  output  is  sampled  "n"  times  more  frequently 
than  the  input  (n  =  2,3,..0)  let  us  adopt  the  notation  that 
Zn  denotes  a  pulse  transform  taken  for  a  sampling  period  of 
T/n.  Then  we  may  defines 


°<Zn>  4  Z;{c<s>}  =£*{»<*>  E<z>} 

H‘zn>  S  {*<=>} 

Once  again  E(z)  is  invariant  (also  see  Reference  4,  page  88) 
so  C(Zn)  ~  H(Zn)*E(z)  but  we  have  a  mixture  of  z  and  Z^ 
transforms,  and  the  meaning  is  unclear,,  If  we  now  define 


then  it  can  be  seen  that  Z  =  zJJ  and  E(Z)  -  E(zjj)  so  in  the 
Zn  "domain"  we  can  write  C(Zn)  =  H(Zn)  E(Z^)C  (2) 

Returning  now  to  the  closed  loop  system  of  Figure  5 
we  can  see  the  following  relationships? 

*<*>’  %8r,  =  ^5?) 

E1(h)  =  R(z)  -  CCz)  (3) 


C  (z)  -  El  (?)  (Zn)‘(?  ^2-n)  j-  (4) 

C  (2*)  =  E1(2)-P(2v.)*  G-CZy,)  (5) 

D/g  ) 

Using  (3)  and  (4):  El(2)  =  - — - - - r  (6) 


Substituting 
(6)  into  (5): 


r(z)  •  2LEz±1  :  (?) 

1  +  &(*■»)} 
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Since  Kdy,)-  there  results:  K(2-»)  =  — —  : 

xt**)  i  +  Zfpt&ycG*)} 


which  can  be  re-written  as: 

KteJ/C {&(*■”)•  £<■*”)} -Ote»).  £•(*■»)  4  K(2»)=0 

Equation  (3)  is  a  functional  equation  of  the  form 

Z{F(a„)}  -  F(2„)  =  -  K (2-k ) 


(8) 


It  can  be  shown 


/9/ 


that  the  solution  is 


F(*0  =  '  -  D(2v,)'  CrC.Z^) 

so  finally,  the  transfer  function  of  the  multi-rate  con¬ 
troller  is: 


D/2  n  -  — ! _J<CgA_ 


(9) 


The  Zn~  transform  of  the  plant  is  simply  obtained  by 
substituting  Zn  for  Z  and  T/n  for  T  in  G( z).  The  Z  -  trans¬ 
form  of  the  overall  process  is  not  so  simply  obtained.  We 
must  start  by  imposing  steady-state  criteria  and  the  re¬ 
straint  of  physical  realizability.  The  following  develop¬ 
ment  is  due  to  Franklin  and  Ragazzini  (Reference  4). 

1 


For  inputs  whose  Laplace  transform  is  R(s) 


„K 


we 


can  write 


(  \  -  Zyi  )  ({—J-  ^ 


C(2v) 

(l-2^)K  K(2„) 

Thus  the  error  sequence,  R(Zn)  -  C(Zn),  can  be  expressed  as 

uo) 

For  a  stable  process  the  Final  Value  Theorem  states  that 

El  (<*>)  =  ,  (i-C)El(^)  UD 

Z-W  +  l 

In  order  for  equation  (10)  to  meet  the  requirements  of 
equation  (11)  we  must  have 

K(Z*)  =  U2> 


(j-z;1)* 
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where  f(Z  )  is  an  arbitrary  function  of  Z  whose  coeffic- 
n  n 

ients  are  to  be  determined. 

Substituting  (12)  into  (10)  results  ins 


s.  PO*)  - 

E1(2")-  (1-2,-r 

and  for  zero  steady  state  error  the  numerator  of  2l(z) 
must  contain  as  a  factor  the  denominator  of  the  inputs 


(13) 


*<V'  *•«.. 


0-^')K-T(2^  =  Q(2„)  (14) 


or  equivalently:  ciW(^(E^) 


=■  O  -for  WisO,i  ...K 
2-m=  1  U5) 


The  restraints  imposed  by  equation  (15)  are  applied  in  order 
to  determine  the  coefficients  of  the  arbitrary  function,, 
f(Zn).  Once  f(sn)  is  known,  equation  (12)  gives  K(Zi). 

The  only  unknown  now  remaining  in  the  fundamental 
design  equation  for  the  digital  controller „  equation  (9), 
is  K(zj^).  This  is  obtained  from  K(Zn)  by  merely  taking  the 
j  powered  terms  of  K(Zn).  .  ..j  =  n,  2n,  3n,  ... 

In  very  special  circumstances,  K(Z^)  may  be  found  in 
a  simpler  way  by  using  so-called  "single-rate1”  techniques 
described  in  Reference  4,  These  special  circumstances  ares 

a.  The  plant  must  be  open-loop  stable. 

b.  K(Zn)  Is  limited  to  a  finite  settling  time  type 


c. 


(i.e. ,  numerator  polynomial  only). 

The  speed-up  factor,  "n‘!,  in  the  controller  out¬ 
put  must  be  selected  so  that  all  transients  at 
output  sampling  instants  are  over  in  one  period 

of  the  input  sampler,  T.  This  is  equivalent  to 
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requiring  that  "n"  equal  the  number  of  terms  in 
the  numerator  of  the  plant # s  pulse  transfer  func¬ 
tion,  G(z). 


If  these  circumstances  obtain,  and  if  the  design  goal  is 
to  achieve  ripple-free  response  to  a  step  input,  then  for 
a  plant  transfer  function  of 


/4/ 

it  can  be  shown'  that  the  digital  controller  transfer 


function  is  given  by: 


The  main  limitation  of  this  design  procedure  using 
single-rate  techniques  is  that  the  speed-up  factor,  "n", 
must  be  the  same  as  the  number  of  numerator  terms  in  G(z), 
which  may  be  out  of  the  question  for  a  high  order  plant* 

For  the  fifth  order  plant  to  be  discussed  later  a  double¬ 
rate  controller  gave  improved  performance,  but  this  con¬ 
troller  could  not  be  designed  by  the  above  procedure*  It 
is  also  worth  noting  that  the  number  of  terms  in  D(Zn)  can¬ 
not  be  less  than  "n". 

Up  to  this  point  the  discussion  of  response  criteria 
has  been  predicated  upon  the  input  being  a  deterministic, 
noise-free  signal,  but  it  also  may  be  desirable  to  judge 
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a  process  on  the  basis  of  its  noise  rejection  capabilities. 
The  following  remarks  concerning  noise  rejection  are  due  to 
Monroe  (Reference  1). 

In  order  to  consider  an  input  signal  which  is  contam¬ 
inated  by  noise,  certain  simplifying  assumptions  are  neces¬ 
sary.  First,  assume  that  the  system  in  the  absence  of  noise 
would  perform  the  ideal  operation  on  the  input,  ie. ,  the 
systematic  error  is  zero.  Second,  assume  the  signal  and 
noise  are  uncorrelated.  Last,  assume  the  noise  is  "white" 
noise  over  at  least  one  sampling  period.  Then  it  can  be 
shown  that  the  mean  squared  error  of  the  process  is: 


*■/ 


% 

where  is  the  mean  squared  value  of  the  input  noise 
CL  ^  is  a  coefficient  of  the  digital  process,  K(z) 

Kte)  =  <2k 


The  "variance  reduction  factor"  is  defined  as  the 
ratio  of  the  mean  squared  error  of  the  output  function  to 
the  mean  squared  value  of  the  input  noise  and  is  seen  to  be 
the  sum  of  squares  of  the  coefficients  of  the  digital  pro¬ 
cess.  Clearly,  the  smaller  this  factor,  the  better  the 
noise  rejection  capabilities  of  the  digital  program  will  be, 
and  in  general  it  is  desired  that  this  factor  shall  be  less 
than  unity.  In  passing,  note  that  the  minimal  prototype 
responses  to  ramp  and  acceleration  inputs  have  variance 
reduction  factors  of  five  and  nineteen,  respectively. 

The  digital  process  can  be  constrained  to  have  no 
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systematic  error  when  the  input  is  taken  as  a  polynomial 
of  a  given  order  or  lower.  For  example,  if  linear  con¬ 
straints  are  imposed,  then  no  systematic  error  exists  for 
ramp  or  step  inputs.  If,  subject  to  linear  constraints, 
the  variance  reduction  factor  is  minimized,  the  following 
expression  is  obtained  for  the  coefficients  of  a  digital 
process  employing  N  samples  of  past  history  and  predicting 
forward  Od  sampling  intervals. 


a(2N-i __  <k  t£L ~  L  1  // 

A/ ( A/ +  l )  A/  (  A/  +  >  )(  aJ  -  /  ) 


and  the  variance  reduction  factor  is: 


#=  c 


The  requirements  that  K(z)  be  able  to  predict  for¬ 
ward  byodT  arises  due  to  the  delay  inherent  in  the  computing 
process,  D(z).  The  discrete  error  value  at  sampling  time  nT 
equals  the  input  signal  at  time  nT  minus  the  output  result¬ 
ing  from  a  control  pulse  applied  to  the  plant  one  sample 
period  earlier,  ie.,  at  time  (n-l)Te  In  order  to  have  zero 
error,  the  output  must  equal  the  input,  so  K(z)  must  predict 
forward  one  sample  period,  thus  ©<  =  1. 

The  second  formula  above  can  be  solved  for  N,  given 
any  desired  value  of  variance  reduction  factor.  In  the  case 
of  unity  variance  reduction  the  result  shows  that  six  samples 
are  required,  and  the  process  pulse  transfer  function  is: 


— 
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It  was  noted  earlier  that  the  presence  of  a  plant  zero 


on  or  outside  the  unit  circle  in  the  z-plane  prevents  the 
realization  of  minimal  prototype  response.  Similarly,  such 
a  zero  will  prevent  realization  of  the  unity  variance  re¬ 
sponse.  Indeed,  when  a  zero  must  be  accounted  for,  it  be¬ 
comes  quite  laborious  to  determine  any  set  of  coefficients, 
a,  ,  that  will  achieve  noise  reduction.  According  to  Monroe 


(Reference  1)  it  is  possible  to  design  K(z)  with  noise  re¬ 
ducing  properties  for  a  plant  with  a  single  zero  outside 
the  unit  circle  by  the  following  procedures 

let  b  =  the  value  of  the  zero  to  be  accounted  for. 


c<  =  unity  for  prediction  forward  on  one  sample 


period. 


The  pulse  transfer  function  of  the  process  using  N  data 
points  is  y_t 

k(*)=  O' 

For  linear  constraints  on  the  process  (ie.,  zero  systematic 
error  for  ramp  or  step  inputs)  then  the  matrix  equation  be- 


These  values  are 


then  employed  in  the  recursion  formula  shown  in  order  to 
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It  was  noted  earlier  that  the  presence  of  a  plant  zero 
on  or  outside  the  unit  circle  in  the  z-plane  prevents  the 
realization  of  minimal  prototype  response.  Similarly,,  such 
a  zero  will  prevent  realization  of  the  unity  variance  re¬ 
sponse.  Indeed,  when  a  zero  must  be  accounted  for,  it  be¬ 
comes  quite  laborious  to  determine  any  set  of  coefficients, 
a^,  that  will  achieve  noise  reduction.  According  to  Monroe 
(Reference  1)  it  is  possible  to  design  K(z)  with  noise  re¬ 
ducing  properties  for  a  plant  with  a  single  zero  outside 
the  unit  circle  by  the  following  procedures 

let  b  =  the  value  of  the  zero  to  be  accounted  for, 
oC  =  unity  for  prediction  forward  on  one  sample 
period. 

The  pulse  transfer  function  of  the  process  using  N  data 

points  is  ... 

/V  K  ,  ,  *p'  -K 

KC  z)= 

xT,  «=' 

For  linear  constraints  on  the  process  (ie.,  zero  systematic 
error  for  ramp  or  step  inputs)  then  the  matrix  equation  be¬ 
low  can  be  solved  for^  3  and/^  .  These  values  are 

then  employed  in  the  recursion  formula  shown  in  order  to 
find  and  so  on  up  to  ^  ^  . 
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The  recursion  formula  is: 


0 
f  M 


H+4 


(/-b+b1)* 

L>  'Mi-7. 


Sven  with  a  noisy  input  signal  it  is  still  desirable 
to  exercise  some  control  over  the  transient  behaviour  of 
the  system.  It  is  possible,  again  according  to  Monroe,  to 
design  the  over-all  process  according  to  a  composite  cri¬ 
terion  wherein  preference  can  be  expressed  for  step  tran¬ 
sient  response  (in  the  mean  squared  error  sense)  over  noise 
rejection  by  means  of  an  arbitrary  parameter, ft  ,  ranging 
between  zero  and  one.  Subject  to  linear  constraints  and 
for  no  plant  zero  outside  the  unit  circle,  the  process  co¬ 
efficients,  a^,  are  calculated  from  the  following  formulae: 
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(^recursion  Ae  /cl±  /  ‘oyi  s h  i p  ) 


0-ux  =  [*'-f-£:/i)ldul  ~ 


d0  =  u'-~ 


&1  = 


Aa  ~  A  / 

/*  4  -  A  4* 


C  $~ta.rt 
’Pormu  I  cl  e) 


A  =  Ahi 

\,-i 


n  =  (N-,)XT'~  NX?  f  X 


A/ 


 (N-/)  X^-  A/  A  7  ■<L  X/ 


A/ 


J. 

2. 


^x0-/jx  (x.-')^ 

=  i.  f 2  tAJ  +  f  AA>12/4  J 

X[=  z(z  }  -  i  {«&)}*[< ^Hth) ] i 


Figure  6  shows  the  Fortran  coding  of  these  equations 
and  the  results  for  a  ten  point  process  with  beta  of  0.75. 
The  program  may  be  used  for  another  design  merely  by  chang¬ 
ing  the  two  statements  giving  the  values  of  the  parameters 
beta  and  AN. 

For  the  design  of  a  process  to  a  composite  criterion 
when  the  plant  has  a  zero  outside  the  unit  circle  see  Ref¬ 
erence  1. 


26 


JOB  STAPP3 

PROGRAM  STAPP3  _ _ 
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Figure  6.  Program  for  computing  composite  coefficients 
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5.  Program  requirements. 

Before  discussing  the  realization  of  the  digital  com¬ 
pensator,  D(z),  by  programming  of  a  digital  computer,  let 
it  be  noted  that  the  pulse  transfer  function  of  the  compen¬ 
sator  can  be  realized  in  at  least  two  other  ways  if  a  digital 

/2/ 

computer  is  not  available  or  economically  feasible.  First 

is  the  realization  by  a  delay-line  network,  sometimes  refer¬ 
red  to  as  a  pulsed-data  processing  unit,  through  the  use  of 
delay  elements,  potentiometers,  and  summing  amplifiers.  The 
pulse  transfer  function  can  also  be  realized  by  a  pulsed- 
data  network  consisting  of  resistors  and  capacitors.  In¬ 
ductors  are  avoided,  because  in  the  usual  frequency  range  of 
digital  and  sampled-data  control  systems  the  sizes  and  weights 
of  the  required  inductors  would  be  prohibitive  and  the  losses 
would  be  excessive.  Pulsed  data  RC  networks  have  the  advan¬ 
tages  of  simplicity  in  structure  and  economy  in  implementation. 

When  a  digital  computer  is  available,  realization  is 
relatively  simple  by  programming.  The  program  can  generally 
be  carried  out  by  means  of  one  of  three  different  methods: 

a.  Direct  programming 

b.  Iterative,  or  cascade,  programming 

c.  Parallel  programming 

Let  the  pulse  transfer  function  (or  Z-transform)  of 
the  digital  controller  be 
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The  first  terra  in  the  denominator  of  D(z)  must  be  a  non¬ 
zero  constant  for  the  sake  of  physical  realizability.  In 
the  direct  method  of  programming  it  is  necessary  that  this 
term  be  unity,  since  the  direct  method  employs  the  differ¬ 
ence  equation  equivalent  of  the  above  pulse  transfer  func¬ 
tion.  Writing  the  difference  equation  we  get 
N-l  M-l 

EoGT  [’nTl  -  Em  ~  2i  Ecur 

j-t 

In  the  digital  computer  this  method  requires  2M+2N-1 
cells  of  storage  for  data  and  constants,  M+N-l  multiplica¬ 
tions,  M+N-2  additions  and/or  subtractions,  and  M+N-2  data 
transfers. 

The  cascade  method  is  based  on  the  fact  that  D(z)  is 
factorable  into  its  poles  and  zeros  such  that  a  cascade  of 
first  order  transfer  function  can  be  formed.  In  the  computer 
this  method  requires  2M+N+3  cells  of  total  storage,  M+N+l 
multiplications,  M+N  additions  and/or  subtractions,  and 
M+2  data  transfers. 

The  parallel  method  of  programming  is  based  on  the 
partial  fraction  expansion  of  D(z)  and  requires  3M+2  cells 
of  storage,  2M  multiplications,  2M-1  additions  and/or  sub¬ 
tractions,  and  M+2  data  transf ers.^^ 

Before  a  choice  of  methods  can  be  made  on  the  basis 
of  storage  requirements  or  computation  time,  D(z)  must  be 
known,  since  any  advantage  of  one  method  over  the  others 
depends  primarily  on  the  number  of  terms  in  the  numerator 
of  D(z).  Generally,  other  considerations  prevail.  For 
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example,  in  the  cascade  method  each  stored  constant  is  either 
a  pole  or  a  zero  of  D(z),  and  thus  this  method  lends  itself 
readily  to  experimental  manipulation  of  the  pulse  transfer 
function  in  a  "cut  and  try"  sense.  If  storage  space  is  at 
a  premium,  as  in  a  special  purpose  digital  computer,  the 
direct  method  enjoys  the  advantage  that  data  transfers  can 
be  effected  by  a  single  command  that  processes  all  the  data 
in  a  shift  register . In  the  simulation  program  described 
in  Appendix  1  the  direct  method  of  programming  the  digital 
controller  is  used. 
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6.  Implementation  of  a  double-rate  digital  controller  by 
direct  programming  methods. 

Direct  programming  takes  the  difference  equation 
equivalent  of  the  pulse  transfer  function  and  interprets 
z"1  to  be  a  delay  of  one  computation  period.  Thus,  in  the 
single-rate  controller  both  input  and  output  discrete  se¬ 
quences  are  equally  spaced  in  time,  and  the  programming  is 
straight-forward.  The  computer  program  for  a  double-rate 
controller  will  be  described  in  terms  of  the  various  arith¬ 
metical  and  shifting  operations  that  must  take  place  and 
not  in  terms  of  any  particular  programming  language.  The 
method  will  be  shown  using  an  example,  and  can  be  easily 
extended  to  the  case  of  an  n-rate  controller. 

Suppose  that  the  pulse  transfer  function  of  the  double¬ 
rate  controller  is: 

0(2Z)  -  4  AT*'?*  +  ^3  +•  +  A  S' •  +  AA  1 

Z~  1+  +- +  6 +  Bte'Zl'*' 

The  output  from  the  controller  will  be  a  sequence,  SoutCZ^.), 
of  discrete  values  separated  in  time  by  T/2  seconds.  The 
input  sequence  is  spaced  at  intervals  of  T  as  determined  by 
the  error  sampling  rate.  When  we  write  the  difference  equa¬ 
tion  to  find  a  new  value  for  Eout  by  cross-multiplying  the 
transfer  function  and  treating  as  a  ^elay  °f  T/2  seconds 

we  get: 

eout  =  Al- +  AZ  •  Ein 

+  A 4  1  £  'Is1  +  &£•  Eim  Eiki 
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SZ*Eout 

+  BS“‘  Eout[(’n-4)t/ji’] 


4-  83  •  Eoot  f  ('N-i)1/aJ  4-  64  Eout  [(M-3)3i] 

4-  Q  b  I  E QOT  ^  -  S’)  J)U~\  I 


Since  the  discrete  sequence,  Ein(z),  is  at  intervals  of  T 
seconds  only  half  of  the  terms  involving  Ein  above  will  be 
non-zero.  Which  half  depends  on  whether  N  is  even  or  odd. 
When  N  is  even  (ie.,  an  output  value  is  generated  at  the 
"same"  time  the  error  is  sampled)  then  the  terms  involving 
Al,  A3,  and  A5  are  non-zero.  When  N  is  odd  (ie. ,  an  output 
value  is  produced  by  the  controller  between  input  samples) 
then  the  only  terms  with  non-zero  values  of  Ein  are  those 
involving  A2,  A4 ,  and  A6.  This  amounts  to  alternately 
multiplying  the  input  sequence  by  the  even  numbered  coef¬ 
ficients  of  the  numerator  of  D^i^),  then  the  odd  numbered 
ones,  ...etc. 

In  a  given  period  of  time  the  output  sequence  from  a 
double-rate  digital  controller  will  consist  of  twice  as 
many  discrete  values  as  the  input  sequence.  For  the  sake 
of  illustration  let  us  number  both  sequences  backward  in 
time  in  the  same  way,  thus:  The  most  recent  value  is  num¬ 
ber  one,  the  next  most  recent  value  is  number  two,  etc. 

The  total  number  of  values  in  each  sequence  that  must  be 
stored  in  the  computer  memory  is  determined  by  the  number 
of  terms  in  the  controller  transfer  function.  In  the  ex¬ 
ample  of  this  section  we  require  output  values  numbers  one 
through  six  and  input  values  one,  two  and  three.  Each  out¬ 
put  value  is  employed  only  once  in  each  numbered  position 
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and  passes  from  storage  after  six  intervals  of  T/2.  Each 
input  value  is  employed  twice  in  each  numbered  position ^ 
once  as  a  multiplier  of  Aj  and  T/2  seconds  later  as  a  mul¬ 
tiplier  of  A.,,.  Since  the  numerator  of  D(Z_)  contains 

J  TJL  £ 

six  terms,  each  input  value  also  passes  from  storage  after 
six  intervals  of  T/2.  A  simple-minded  pictorial  represen¬ 
tation  may  serve  to  clarify  these  ideas.  Let  the  present 
time  be  the  16th  instant  of  the  closing  of  the  output  samp¬ 
ling  switch  of  the  double-rate  controller,,  Neglecting 
computation  delay  (which  is  usually  at  least  an  order  of 
magnitude  less  than  the  sampling  period)  the  present  time 
is  also  the  eighth  closing  of  the  input  sampling  switch. 

The  input  and  output  sequences  may  be  thought  of  as  exist¬ 
ing  in  storage  as  shown  in  Figure  7. 

Input  sequence,  Ein 


observation 

©^ 


© 
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8T 


present 

time 


7T 


CT 


time 


Output  sequence,  Sout 

®  © 


© 


to  be  computed 

/  t  t 

16T/2  15T/2  14 T/2  13T/2  12T/2  11T/2 


time 


(older  values) 

Figure  7,  Input  and  output  sequences  for  EK^) 


In  the  difference  equation  the  terms  involving  the  input 
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sequence  would  be  (since  in  this  example,  N  =  16): 
Al*Ein(8T)  or  Al*Ein(l) 

A3*Ein(7T)  or  A3*Ein(2) 

A5‘Ein(6T)  or  A5*Ein(3) 

The  value  of  Ein(l)  has  just  been  measured  by  the  input 
sampler  and  the  controller  produces  Eout(l). 

Now  let  us  move  on  in  time  by  T/2  seconds.  It  is 
necessary  to  re-number  the  output  sequence  since  a  new 
Eout(l)  will  be  produced.  The  previous  Eout(l)  becomes  the 
current  Eout(2),  etc.,  with  the  previous  Eout(6)  being  dis¬ 
carded.  We  are  now  between  input  samples  and  no  new  input 
information  is  to  be  received.  The  previous  input  values 
have  become  T/2  seconds  older,  however,  and  the  sequences 
in  storage  might  be  pictured  as  shown  in  Figure  8e 

Input  sequence,  Ein 


(D 


8T 

to  be  computed 

d/  ® 


17 T/2  IS T/2 

present 

time 


7T 

Output  sequence,  Eout 

©  0 


6T  time 


15T/2  14T/2  "T3T/2  “  T2T/2 

(older  values) 


time 


Figure  8.  Sequences  between  input  samples 


Since  the  difference  equation  is  in  terms  of  T/2  the  non-zero 
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terms  involving  the  input  sequence  have  changed  and  are 
now  (for  N  =  17): 

A2»Ein(8T)  or  A2»3in(l) 

A4*  Ein ( 7T)  or  A4*Sin(2) 

A6*Ein(6T)  or  A6*Ein(3) 

The  controller  computes  a  new  output  value  using  all 
of  the  denominator  coefficients  of  its  transfer  function 
but  only  the  even-numbered  numerator  coefficients.  As 
before,  the  output  sequence  must  be  shifted  since  we  are 
only  half  a  period  away  from  receiving  a  new  observed  value 
for  Ein(l). 

For  an  example  of  the  implementation  of  these  ideas 
using  the  CDC  Fortran  programming  language  see  Appendix  1. 


35 


7.  Applications. 

The  results  to  be  described  in  this  section  were 
obtained  using  the  digital  computer  program  simulation  of 
a  control  system  described  in  Appendix  1.  The  response 
curves  were  obtained  through  use  of  a  Graphplot  computer 
routine  for  the  CDC  160  digital  computer  devised  by  Lt. 
Robert  L.  Hogg  of  the  U.  S,  Naval  Postgraduate  School , 
Monterey,  California.  The  program  in  Appendix  1  also  pro¬ 
duces  numerical  results  on  a  printer,  but  these  are  not  in¬ 
cluded  in  this  work. 

Initial  investigation  was  made  of  a  simple  plant 
described  by  a  second  order  differential  equation  whose 
pulse  transfer  function  had  no  zero  outside  the  unit  circle 
in  the  z-plane.  Next,  a  third  order  plant  with  a  pulse 
transfer  function  containing  a  zero  outside  the  unit  cir¬ 
cle  was  investigated.  Both  of  these  plants  were  theoretical 
and  in  neither  case  was  an  actual  physical  counter-part  at 
hand. 

At  the  U.  S.  Navy  Electronics  Laboratory,  San  Diego, 
California,  a  group  headed  by  John  B.  Slaughter  is  engaged 
in  studies  of  Digital  Sampled-Data  Feedback  Systems  for 
Shipboard  Control  Applications  (NEL  problem  N4-17).  The 
author  was  privileged  to  spend  ten  weeks  with  that  group 
early  in  1963,  during  which  time  a  twin  3"/50  gun  mount 
was  in  the  process  of  being  installed  to  serve  as  the  plant 
in  an  experimental  control  system.  A  USQ-20  general  purpose 
digital  computer  is  to  be  employed  as  the  controller.  The 
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manufacturers  of  the  gun  mount  supplied  a  fifth  order 
transfer  function  for  the  gun  in  train.  The  majority  of 
the  simulation  studies  to  be  described  in  this  thesis  were 
conducted  using  this  transfer  function,  and  it  is  felt  that 
the  results  obtained  are  applicable  not  only  to  this  par¬ 
ticular  obsolescent  plant  but  also  to  present  day  ship¬ 
board  weapons,  especially  missile  launchers. 

The  simulation  program  has  the  capability  of  testing 
a  system  for  input  signals  of  either  unit  step,  ramp,  or 
acceleration.  In  addition,  additive  Gaussian  noise  of 
selected  mean  and  variance  may  be  included  m  the  input. 

The  first  action  of  the  program  is  to  produce  the  response 
obtained  in  the  case  of  the  continuous,  un-compensated  sys¬ 
tem  for  purposes  of  comparison.  In  the  event  this  response 
is  unstable  no  graph  is  made  of  it.  Instead,  the  test  input 
is  reproduced.  The  next  action  of  the  program  is  to  simulate 
the  same  plant  incorporated  in  a  sampled-data  system  with  a 
single-rate  digital  controller  and  produce  the  response  for 
this  case.  The  program  is  further  capable  of  simulating  a 
system  using  a  double-rate  digital  controller  to  produce  a 
third  response  curve,  although  this  is  not  done  in  every  ex¬ 
periment.  For  convenience,  the  program  will  also  compute 
the  pulse  transfer  function  of  the  controller  when  requested 
to  do  so,  or  the  controller  coefficients  may  be  provided  by 
the  user.  The  Z-transforms  employed  in  this  work  were  com¬ 
puted  by  the  program  described  in  Appendix  2  and  were  veri¬ 
fied  by  hand  for  the  second  and  third  order  systems. 
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The  second  order  plant  has  a  Laplace  transfer  func¬ 
tion  given  by  i  / /^\  _  //  41 

HCS  ~  S(s-hZS) 

The  pulse-transfer  function  of  the  plant  and  zero-order  hold 
combination  for  a  sampling  period  of  T  =  .01  seconds  is 

r/z)  -  .0£ajZ£(l+  d. 

(i--2~')(/-o-~7i 2*oi  a~‘)  ~  /- /. 77220/ h'V  o.nzqoi  z~‘ 

There  are  no  zeros  of  G(z)  on  or  outside  the  unit  circle  in 

the  z-plane,  so  minimal  prototype  compensation  can  be  a~ 

chieved.  As  noted  earlier ,  the  minimal  prototype  over-all 

pulse  transfer  function  for  a  step  input  is  K(z)  =  z“^ ,  so 

I  __  /-  /,-77  ggo/e‘V  .JJJJCjl  g~' 

'  <?(2)  I  -  « (2)  ~  .  as’222~‘  +-  *64  86  2  2~‘ 

2 -  /£73f3fry-/#  7^0/92"' 
•  Og'LB^'Cl  +  ,cn.G>4SSZ~')  /.o  -h  0.9204  S’-TZ"' 

(note  that  l-K(z)  contains  the  pole  of  G(z)  lying  on  the 
unit  circle. ) 

Figure  9  shows  the  plant  behaviour  in  a  continuous,  un¬ 
compensated  system  and  the  behaviour  in  a  sampled-data 
system  using  the  digital  controller  above.  Note  that  there 
is  zero  systematic  error  at  the  first  sampling  instant  after 
the  input  was  applied,  and  also  at  all  succeeding  sample 
instants.  For  a  K(z)  containing  "n"  terms  there  will  be 
zero  error  for  a  design  input  after  "n"  sampling  instants. 

The  plant  behaviour  between  samples  is  another  matter,  and 
in  this  example  probably  would  not  be  considered  satisfac¬ 
tory,  unless  the  only  desire  was  to  increase  the  speed  of 
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response  with  no  regard  for  over-shoot  and  inter-sample 
ripple. 


Simulated  Gaussian  noise  with  a  standard  deviation 
of  0.204  was  added  to  the  input  and  the  experiment  was  re¬ 
peated.  Figure  10(a)  shows  the  results.  According  to  Mon¬ 
roe  the  process  should  have  unity  variance  reduction  factor 
since  the  sum  of  the  squares  of  the  process  coefficients  is 
one.  Judging  from  Figure  10(b),  which  shows  the  noisy  test 
input,  this  is  approximately  correct,  since  no  noise  reduc¬ 
tion  is  apparent. 

In  an  effort  to  improve  the  transient  behaviour  for  a 
step  input  a  staleness  factor  of  0.5  was  selected.  The 


over-all  process  was  designed  from 


K(*)  - 


a-  YU  rrtt  'yntL  / 


/-  O.S  2"' 

and  the  Final  Value  Theorem  showed  that  for  zero  steady- 


state  error,  a  =  l-  c  =  C.5  so  that  K(s:)  = 


O. 


-/ 


/  -  O.  S’. 


D(2)  = 


,  _  OS-Z-' 

(/-■Z  'J(/-0.77S*0/S  J  — 


.  OS3X  E~'+  .  04 2^  S' 


/ 


/-  O.J-Z 


_  0.32CHOOS’2 

S-  .  04202.-' 

The  step  response  using  this  digital  controller  is  shown 
in  Figure  11.  This  is  an  example  of  a  non-finite  settling 
time  process,  but  for  practical  purposes  zero  error  is 
achieved  at  about  the  seventh  sampling  instant  and  there¬ 
after,  and  the  inter-sample  ripple  and  over-shoot  is  much 


<2  40W7  —  7 374‘aa 
/*Q  -f-  0 .9 304 Z~' 
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reduced  compared  to  the  minimal  design. 

This  same  staleness  factor  design  was  subjected  to 
a  noisy  step  input  with  the  results  shown  in  Figure  12. 

Since  K(z)  can  be  expressed  as  an  infinite  series  by  long 
division,  thus:  0.5z-1  +  . 25z~^  +  .125z~3  +  . 0625z~4. . . . 
the  variance  reduction  factor  is  approximately  one-third. 

It  is  not  apparent  from  the  simulation  results  that  the 
process  indeed  had  this  noise  reducing  property. 

The  six  point  finite  memory  process  for  unity  variance 
reduction  (see  page  23  above)  was  then  subjected  to  a  step 
input  with  and  without  noise.  Figures  13  and  14  show  the 
results.  This  process  and  the  following  one  depend  on  the 
pulse  transfer  function  of  the  plant  not  having  a  zero  out¬ 
side  the  unit  circle.  It  can  be  seen  in  Figure  13  that  the 
transient  behaviour  of  a  design  whose  sole  criterion  is 
variance  reduction  is  rather  poor.  For  this  reason  the  next 
process  considered  v/as  a  composite  one  in  which  three-quarters 
of  the  design  effort  was  towards  transient  control  and  only 
one-quarter  towards  noise  reduction.  The  program  on  page 
27  was  used  to  compute  the  process  coefficients,  as  shown 
on  the  same  page.  This  gave  the  most  complex  controller  so 
far  with  11  numerator  terms  and  11  denominator  terms.  The 
results  of  the  simulation  experiment  with  this  controller 
are  shown  in  Figures  15  and  16.  The  transient  behaviour 
displayed  in  Figure  15  is  a  noticeable  improvement  over  that 
of  Figure  13  where  the  only  aim  was  noise  reduction. 

At  this  point  a  certain  disenchantment  with  the  notion 
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of  noise  reduction  begins  to  set  in.  The  plant  is  a  simple,, 
well-behaved  one,  and  in  the  four  experiments  so  far  using 
"clean’1  input  signals  both  the  continuous  system  responses 
and  the  sampled-data  system  responses  can  be  verified  ana- 
lytically  with  some  effort.  When  the  unit  step  input  has 
"white"  noise  added  to  it  as  shown  in  Figure  10(b)  two  ten¬ 
tative  conclusions  suggest  themselves.  First,  the  contin¬ 
uous  system  with  no  compensation  is  relatively  unaffected 
by  the  presence  of  the  noise.  A  close  inspection  of  Figure 
10(a)  reveals  that  the  peak  overshoot  still  occurs  at  about 
one-tenth  of  a  second  and  its  value  is  increased  only  slight¬ 
ly  (from  29%  to  37%).  A  steady  state  error  equal  to  the  mean 
squared  value  of  the  noise  can  be  inferred  from  Figures  12 
and  16  in  which  the  time  scale  has  been  doubled  for  a  longer 
look  at  the  plant0 s  response.  Over-all ,  however,  the  low- 
pass  character  of  the  plant  effectively  prevents  any  drastic 
reaction  to  the  presence  of  the  noise.  Not  so  in  the  case  of 
the  sampled-data  system,  so  our  second  thought  at  this  point 
is  that  the  insertion  of  the  digital  controller  in  the  for¬ 
ward  path  disrupts  the  low-pass  filtering  action  of  the  plant. 
The  noise  in  the  error  signal  is  the  same  as  in  the  input,  and 
the  noise  power  in  each  discrete  error  sample  is  integrated 
over  the  next  sampling  period  and  applied  to  the  plant.  This 
integrated  noise  power  causes  some  fairly  drastic  gyrations 
in  the  plant  regardless  of  the  design  criteria.  A  brief  com¬ 
parison  of  the  four  noisy  experiments  up  to  this  point  provides 
small  reason  for  employing  a  variance  reduction  design  in 
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preference  to  a  noise-free  criterion.  It  appears  that  the 
effort  could  better  be  expended  in  pre-filtering  or  other 
measures  to  eliminate  noise  prior  to  the  input  to  a  sampled- 
data  system. 

The  next  four  simulation  studies  were  conducted  using 
a  third  order  plant  whose  transfer  function  is  given  by 


- - 

s(  s  s-  as-J(s+  70) 


For  T  =  .01  seconds  the  z-transform  of  the  plant-and-hold  is 


.  O/O  &O0-2  .  0336  72.  -t-  .OQ6~S~g^/S  3 


/.O  -  3.  37S32C>-Z~‘+  /  6i^72"X-^3?67^/2  3 
or  in  factored  form 

a.  967  376~  2~' ) 

(/-2~)  (/  -  0.7  7226/Z~'X/-  B~') 

The  presence  of  a  zero  of  G(z)  outside  the  unit  circle  is 
noted.  This  zero  will  prevent  us  from  realizing  a  "minimal" 
prototype  response.  In  accordance  with  the  rules  given  on 
page  16  the  design  of  a  process  giving  prototype  response 
to  a  step  input  proceeds  as  follows: 

KU)  -  (/  +  3'90>727<s-  2~') 

!-/<(*)  —  (/-  2~')(/+  ^  2~'J 

Solving  for  the  two  unknowns,  and  b^#  results  in 
CL±*  .2&OL7.  bt^.14193%  and  thus 

K(  2)  =  O.  31 S  20  6  2  2_/  ~h  O.  7-/7  93  5?S-=l 

=  ft  +  3.  16  7  27J-Z~')0‘  33-30622-' 

Solving  for  the  digital  filter  coefficients  by  hand,  we  write 
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tf>Z^(/-,778$0/B')(/-.496S?S'2')  JS20&2^' 

^  =  ('.0/at)(/+.  309327j/Xte^J¥*:Kr7Z^Cr^^ 

and  after  cancellations  of  the  factors  shown,  we  get 

9(2)  =  23.719934  -  3Q.3^7f?3£~y 

/.O  i-  0,  9S73&-7 E~'  9  0.  /sT^  2  2“  a 

The  un-compensated  response  and  the  response  of  the  sampled- 
data  system  using  the  D(z)  above  are  shown  in  Figure  17. 

The  non-minimal  nature  of  the  digitally-controlled  response 
is  characterized  by  the  fact  that  it  took  two  sampling  per¬ 
iods  (.02  seconds)  for  the  output  to  reach  zero  error  at 
sampling  instants.  There  is  a  slight  amount  of  inter¬ 
sample  ripple,  but  the  peak  overshoot  is  greatly  reduced 
from  what  it  was  for  the  continuous  system.  The  improved 
damping  of  this  third  order  sampled-data  system  compared  to 
the  second  order  system  considered  previously  is  attributed 
to  additional  filtering  action  of  the  pole  at  -70.  Since 
this  design  was  for  a  step  input  there  will  not  be  zero  sys¬ 
tematic  error  for  higher  order  inputs,  such  as  a  ramp.  This 
can  be  observed  in  Figure  18  where  the  curve  labelled  D1  is 
the  response  of  the  digitally  controlled  system  to  a  ramp 
input.  (Note:  in  this  figure,  and  all  others  for  which  the 
test  input  is  a  ramp,  a  45  degree  sLope  from  the  origin  is 
the  input  signal).  The  controlled  response  settles  down  to 
a  steady  state  error  of  about  .017  units  after  two  sample 
periods.  Whether  or  not  this  is  an  intolerable  error  in  re¬ 
sponse  to  a  velocity  input  depends  on  the  system.  For  ship¬ 
board  fire-control  application  such  as  tracking  a  target,  a 
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design  for  zero  error  to  a  ramp  is  usually  preferable. 

More  about  this  later. 

An  additional  refinement  of  the  step  response  of  this 
third  order  system  can  be  achieved  at  the  cost  of  one  extra 
period,  T,  in  reaching  steady  state  by  designing  for  ripple- 
free  response  as  discussed  on  page  17.  For  this  case  we  get 

tffr)  =  (/  r  O-  3093  27  S~')(/  ■+ 

-  (/-  *~'X/  y-  4'5'/  +  ^a.  Z~Z) 

and  by  equating  coefficients  of  like  powers  of  z'k  (k=l,2,3) 
we  obtain  the  below  equations  to  solve  for  al,  bl,  and  b2. 

-  /  -  b± 

^  •  /  ~7  Cs  Q  2.  ^2^  ~  ^ 

0.  62//30  2  O.JL  = 

Solving  these,  we  get  al  =  .2084318  so  that  the  process  is 
K(z)  =  . 2084318z-1  +  ,662105z“2  +  . 12946 34z~3 

These  process  coefficients  as  well  as  the  pulse  trans¬ 
fer  function  of  the  plant  were  provided  to  the  simulator 
program  which  then  calculated  the  digital  controller  and 
produced  the  ripple-free  response  shown  in  Figure  19. 

Figure  20  shows  the  same  system  with  a  velocity  input.  The 
steady  state  ramp  error  has  increased  slightly,  from  .017  to 
.020  due  to  the  additional  term  in  the  process  z-transform, 
and  the  controller  had  a  total  of  six  terms.  In  the  absence 
of  input  noise  we  have  demonstrated  satisfactory  digital  con¬ 
trol  of  this  lightly  damped  plant  when  the  input  is  a  unit 
step. 
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Figure  10(a) 
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Turning  now  to  the  twin  3"/50  gun  mount*,  the  transfer 
function  of  the  forward  path  train  drive  is  given  by 
U  (s )  _  _ /•  60  &  4  a.3  £~  x  /G? _ 

S  +//XS+  /7& S  +  33X3+  33) 

Most  of  this  very  large  forward  gam  is  supplied  by  the 
amplidyne  drive  units  and  results  in  a  quickly  responsive 
system  when  suitable  compensation  is  added.  In  actual  use 
the  train  drive  function  was  elaborately  compensated  by 
various  continuous  networks  in  a  multitude  of  feedback 
paths  in  order  to  obtain  stable  operation.  All  of  this 
continuous  compensation  has  been  removed  from  the  gun  by 
the  group  at  NEL*  and  the  necessary  control  will  be  achieved 
by  digital  methods.  Figures  21  and  22  were  provided  by 
Donald  Lackowski*  an  engineer  in  that  NEL  group*  and  show 
the  root  locus  and  Routh  criterion  analysis  for  the  forward 
transfer  function  with  unity  feedback.  It  cart  be  seen  that 
the  stability  limit  on  the  gain*  K*  is  roughly  three  orders 
of  magnitude  lower  than  the  actual  value  present*  and  thus 
the  plant  is  highly  unstable  without  compensation.  This  is 
borne  out  by  the  simulation  studies  to  be  described  next. 

In  every  case  the  continuous  response  went  rapidly  off  scale 
and  could  not  be  shown  in  the  figures  to  follow.  As  previ= 
ously  mentioned*  in  this  event  the  simulator  program  re 
produced  the  test  input  in  lieu  of  the  unstable  response. 

The  computer  program  described  in  Appendix  2  was  ex¬ 
ceedingly  helpful  at  this  point  since  finding  the  z-trans- 
form  of  this  fifth  order  equation  for  just  a  single  value 
of  T  is  quite  laborious  by  hand.  With  the  computer  a 
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number  of  transforms  was  quickly  obtained  for  a  spectrum 
of  sampling  rates  ranging  from  one  hundred  samples  per 
second  to  as  slow  as  one  sample  every  four  seconds*  For 
the  faster  rates  the  pulse  transfer  functions  themselves 
were  of  the  fifth  order  (ie.„  five  numerator  terms  and  six 
denominator  terms).  As  the  rate  was  decreased  it  was  ob¬ 
served  that  the  higher  ordered  terms  in  both  numerator  and 
denominator  of  G(z)  became  progressively  smaller*  For  T 
of  one-half  second  or  larger  a  good  approximation  of  G(z) 
can  be  made  using  only  a  second  order  function*  In  addi¬ 
tion,  the  zero  of  G(z)  which  lay  outside  the  unit  circle 
for  smaller  values  of  T  had  moved  inside  the  circle  for 
T  =  0.5  seconds.  The  computer  calculated  pulse  transfer 
function  of  the  plant  and  hold  combination  for  T  of  one- 
half  second  is  given  below: 


Numerator  of  G(z) 

Denominator  of 

1336.33153319Z"1 

1.00000000 

897. 96714053z“2 

-1. 004300 37 z"1 

6.02702748z~3 

*00430125z"2 

. 00049 220z"4 

-.0000Q087z~3 

.00000000z“5 

.  00000000 z**4 
.00G000Q0z“5 

The  computer  program  also  factored  the  numerator  and  gave 
the  following  real  values  for  the  zeros  of  G(z)? 
-.00000001,  -.00008267,  -.00669675,  -.66518491,  0*0 
It  was  decided  to  investigate  the  advantages  of  a  multi- 
rate  sampled  data  system,  the  analysis  and  design  of  which 
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was  discussed  beginning  on  page  17.  The  simulator  pro¬ 
gram  has  the  ability  not  only  to  produce  results  for  the 
continuous  system  and  a  single-rate  digitally  controlled 
system,  but  also,  upon  request,  to  simulate  a  system  using 
a  digital  controller  whose  output  to  the  plant  occurs 
twice  as  frequently  as  the  input  error  samples  are  sup¬ 
plied  to  the  controller.  This  is  called  a  double-rate 
controller.  The  output  rate  was  chosen  to  be  two  samples 
per  second  in  order  to  take  advantage  of  the  second  order 
approximation  to  G(z)  mentioned  above.  It  was  desired 
that  the  over-all  process  have  zero  steady  state  error  to 
a  velocity  (ramp)  input.  Following  the  method  of  page  19 
the  design  for  prototype  ramp  response  proceeds  as  follows s 
The  fundamental  design  equation  for  the  double-rate  con¬ 


troller  is 


For  T/2  =  0.5  seconds,  G(Zy)  ~  g— 

^  /  —•  if  ^2 


Applying  steady  state  and  realizability  constraints,  we  get 


and  a  look  ahead  at  the  constraint  equation, 


Hj'l 
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shows  the  need  for  only  two  free  constants,  a,  and  a-*  so 

0(2*)  =  "  Pfei.  )*  F(^x)  -T*£a  (Q±2z' + 

=  o.rzj'-z^jzr+^A  =  o.&;‘- a^-*.  a3z\ 


■for  m  =  0  :  QCH^) 


2,  =  1 


-  0,S 


C !j_  —  <0^  —  o 


■for  ~m  =  1  t  dQ 

a?;' 


=  o.s-  -  ?a±-  4  a*.  =  o 


Za*l 


solving  these  equations  we  gets 


1.5,  a  — 


-1.0 


so  that  K(za)=  (i^js )  (l.rz-'-2r)=  (i  +  z^i-rej'-z;1 

=  A<3 

To  find  K(2z)  we  take  the  even -powered  terms  with  the  re- 

suit  K(z*)=  3.  o  -  /.  o  z~  *  =  /cY^J 

which  can  be  recognized  as  the  minimal  prototype  ramp 
response  function  shown  on  page  14  expressed  in  the 
domain. 

Substituting  now  in  the  fundamental  design  equation, 
the  pulse  transfer  function  of  the  double-rate  controller 

n,  x  (/'  2*)  (/•  *  2-°  ** a-a^ga3  - 

^  2  /  (/534>b;'+  M22*)  (  /  -  3.0  -h  /.0  ?/■  *) 

=  /  ^-  -  gj7 _ 

/336  -  43  2  Z^  -  2“* 

We  must  normalize  this  result  with  respect  to  the  leading 
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denominator  coefficient,  since  direct  programming  requires 
a  value  of  unity  in  this  position.  The  final  result  is 


DteJ  = 


.60//a27S-4  -  .  0007  4  $  S  O 
/.0  -  0.327 2443!  —  0.672/^69  Z "  x 


Before  proceeding  with  the  simulation  run  we  must 
provide  for  the  single-rate  action  with  T  =  1.0  seconds. 
This  is  simply  done  by  providing  the  simulator  program  with 
the  ramp  minimal  prototype  function  and  G(z)  coefficients. 
For  this  sampling  rate  an  even  better  second  order  approx¬ 
imation  is  available  for  G(z),  namely? 

,  ,  ,  _  3S-7&  3~'  +  924  g~a 

6rt*J_  ,  _  _-/ 


The  responses  of  the  two  systems  for  a  noise-free 
velocity  input  can  be  seen  in  Figure  23.  The  single-rare 
controller  gave  zero  error  at  the  second  sampling  instant 
and  only  moderate  overshoot  and  inter-sample  ripple.  The 
double-rate  system  was  0. 5  seconds  faster  but  over-shoot 
was  worse  due  to  the  larger  control  forces  involved.  The 
ripple  in  the  double-rate  case  was  more  persistent,  also. 

These  same  systems  were  then  subjected  to  step  inputs 
with  the  results  displayed  in  Figure  24.  The  highly  tuned 
nature  of  any  prototype  design  is  exhibited  here,  where 
the  input  is  of  lower  order  than  that  for  which  the  system 
was  designed.  An  interesting  result  in  Figure  24  which 
was  not  expected  is  the  slightly  decreased  peak  over¬ 
shoot  in  the  case  of  the  double-rate  controller. 

The  next  design  criterion  applied  was  ripple-free 
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step  response  using  both  single  and  double-rate  control¬ 
lers  and  the  same  error- sampling  rate  of  the  preceding 
design.  The  ripple-free  process,  K(z),  was  computed  and 
supplied  to  the  simulator,  which  calculated  the  single¬ 
rate  controller,  Dl.  The  double-rate  controller  was  de¬ 
signed  using  the  techniques  discussed  beginning  at  the 
bottom  of  page  20.  Both  controllers  were  quite  simple 
and  are  given  below: 

[)l  -  .  Q2  _  »QQQ4475~--  ,  00000/ 43  d- 

I  +  0-los- 3333  2”'  9  AO 

Figure  25  shows  both  systems  responding  without  ripple  to 
a  noiseless  step  input.  Again,  the  double-rate  system  is 
about  one-half  a  second  faster.  The  ramp  response  is  next 
shown  in  Figure  26.  Both  systems  have  a  steady  state  error 
since  they  were  designed  for  a  step  input  only.  In  the 
case  of  the  double-rate  controller  the  sustained  ripple 
is  due  to  the  fact  that  the  design  technique  was  a  special 
one,  explicitly  for  ripple-free  step  response,  and  cannot 
be  expected  to  give  conventional  behaviour  for  a  different 
input. 

The  results  of  the  simulation  studies  using  the  second 
order  approximation  designs  with  inputs  contaminated  with 
Gaussian  noise  tend  to  reinforce  the  tentative  conclusions 
mentioned  earlier.  It  was  impossible  to  observe  the  response 
of  the  continuous  system  to  a  noisy  input  since  che  contin¬ 
uous  system  was  unstable.  The  disturbing  result  is  that 
the  sample-data  systems,  bcr,h  with  single  and  double-rate 
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digital  controllers,  which  had  given  stable  performance 
of  varying  degrees  of  goodness  when  the  input  was  noise- 
free  were,  with  but  a  single  exception,  unstable  in  the 
presence  of  noise.  The  exception  is  shown  in  Figure  27 
in  which  the  design  of  D1  and  D2  was  for  ripple-free  re¬ 
sponse  to  a  step  and  the  input  was  a  noisy  step.  It  seems 
apparent  at  this  point  that  a  large  signal-to-noise  ratio 
is  a  sine  qua  non  in  the  input  to  a  sampled~data  system. 

In  the  preceding  experiments  it  was  only  the  pulse 
transfer  function  of  the  plant  and  hold  combination  that 
was  approximated  in  order  to  simplify  the  design  calcula¬ 
tions.  The  differential  equation  of  the  plant  itself 
remained  of  the  fifth  order  for  the  numerical  methods 
solution  of  all  response  values. 
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The  simplicity  of  the  design  techniques  for  sampling  rates 
of  one  or  two  samples  per  second  which  result  in  valid 
simplifying  approximations  is  quite  appealing,  but  there 
lurks  an  intuitive  suspicion  that  these  rates  are  too  slow 
for  a  system  required  to  accurately  follow  a  high-speed 
aerial  target.  This  suspicion  has  been  confirmed  in  pri¬ 
vate  communications  between  John  3.  Slaughter  and  Gene  F. 
Franklin,  the  results  of  which  are  that  the  damped  natural 
frequency  of  the  gun  mount  in  train  is  about  one  and  one-half 
cycles  per  second  and  that  sampling  should  occur  between  six 
and  twelve  times  the  damped  natural  frequency  of  any  system 
for  good  results.  For  these  reasons  a  sampling  rata  of  ten 
samples  per  second  has  been  decided  on  by  the  group  at  N3L, 
and  the  remaining  simulation  studies  of  this  work  will  use 
this  rate. 

For  T  -  0.1  seconds  the  pulse  transfer  function  of 
the  plant  and  zero-order  hold  is  shown  below; 


Numerator  of  G(z) 

Denominator  of 

21. 91 241 531 z"1 

1.00000000 

131. 76730855Z-2 

-1. 61605197z-1 

63. S4130401z~3 

. 72871393z~^ 

3.03436539Z-4 

-.11377672Z-3 

.00419043Z-5 

,00612622z-‘i 

-.  00000152z"'J 

The  roots  of  the  numerator  are  located  at  the  origin  and  at. 
-.00139903,  -.05265621,  -.47316625,  and  -5.48614151 

The  presence  of  the  zero  of  G(z)  outside  the  unit 


circle  complicates  the  design  effort  slightly  as  it  did 
in  the  third  order  plant  considered  earlier.  The  same 
methods  apply  and  the  resulting  pulse  transfer  function 
of  the  process  which  gives  prototype  response  to  a  step 
input  is  K(z)  =  .154174866Z"1  +  . 845825134z”2 

Figure  28  shows  this  process  responding  to  a  noise- 
free  step  input,  and  Figure  29  shows  the  same  process  with 
a  noise-free  velocity  input.  The  remarks  made  earlier  con¬ 
cerning  prototype  processes  apply  and  will  not  be  repeated 
here.  Figure  30  shows  the  results  obtained  with  this  de¬ 
sign  when  the  step  input  contains  additive  noise.  The 
best  that  can  be  said  about  the  plant's  behaviour  in  this 
case  is  that  it  remained  stable  and  in  the  vicinity  of  the 
input  signal.  Not  so  for  a  noisy  ramp  input,  to  which  the 
response  of  the  sampled  data  system  was  unstable.  It  ap¬ 
pears  that  the  sensitivity  of  a  sampled  data  system  to 
noise  increases  with  the  order  of  the  input  signal. 

The  prototype  design  for  zero  steady  state  error  in 
response  to  a  velocity  input  begins  with  the  equations  below: 
K(z)  =  (1  +  5.4861415z"1)(a1z"1  +  a2z”2) 

1  -  K(z)  =  (1  -  z_1)2(l  +  bjZ-1) 

From  these  are  obtained  three  equations  to  be  solved  for 
three  unknowns ,  namely 

a1  =  2  -  b1#  -5.4861415a2  =  b^#  and 
a2  -  5. 486141 5ax  =  2h±  -1 

Solving  these,  the  resulting  over-all  process  is 
K(z)  =  . 438754756z_1  +  2. 122490488z"2  -  1. 561245244z~3 
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The  coefficients  of  this  process  as  well  as  the  coeffi¬ 
cients  of  the  numerator  and  denominator  of  G(z)  and  the 
value  of  the  zero  outside  the  unit  circle  were  supplied 
to  the  simulator  program,  which  then  calculated  a  single¬ 
rate  digital  controller.  The  responses  of  the  system  to 
a  noise-free  ramp  and  a  noise-free  step  are  shown  in  Fig¬ 
ures  31  and  32  respectively.  With  the  addition  of  noise, 
instability  again  resulted  when  the  input  was  a  rarnp,  while 
for  a  noisy  step  the  output  was  at  least  bounded  for  the 
duration  of  the  experiment,  although  it  began  to  look,  in¬ 
auspicious  towards  the  end.  See  Figure  33.  In  an  attempt 
to  achieve  at  least  a  semblance  of  stability  in  response 
to  a  noisy  ramp  signal  the  noise  deviation  was  reduced 
from  0.2  to  .05,  but  to  no  avail.  A  future  investigation 
will  be  made  into  the  stability  threshold  for  noisy  ramp 
signals. 

It  was  then  decided  to  design  for  prototype  step  re¬ 
sponse  using  a  staleness  factor  of  0.5.  The  design  steps 
are  shown  below: 

Prototype  K(z)  =  .  154174866Z-'*'  +  .  845825134z“2  as  found 

earlier.  This  function  contains  the  zero  of  G(z)  located 

outside  the  unit  circle  at  z  =  5.48614151.  The  process 

pulse  transfer  function  with  a  0.5  staleness  factor  Is 

given  by  /-z  >>  _  ^2  (o^JSj4jn_4  &  2  ‘jt-QiJLdJL (34j_2_ 

~  /  -  a,  sTH-' 

For  a  unit  step  input,  application  of  the  Final  Value  Theo¬ 
rem  shows  that  "a"  must  have  the  value  0.5.  Thus  K(z)  is; 
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Kfc)  = 


,^77^^33  Z-/  ^ 

/  -  &.  S’  Z~' 


Substituting  in  D/s)  =  ^  (^.). 

^  J  P<z)  I-KC2) 

where  P(z)  is  the  numerator  of  G(z)  given  on  page  71, 


and 


Q(z)  is  the  denominator  of  G(z)  given  on  page  71, 

—  1  1 

and  cancelling  the  common  factors  (1-z  )(l+5.48614151z  ), 

the  single-rate  digital  controller  is  given  by 


Numerator  of  D(z) 

Denominator  of 

. 003517980 

1.00000000 

-. 002167290z_1 

. 95013414z-1 

, 000396307z-2 

. 24861928z-2 

-. 000021510z“3 

.  01088370z"”3 

.000000025z"’4 

. 00001039Z-4 

This  pulse  transfer  function  was  supplied  to  the  simu¬ 
lator  program?  test  inputs  of  noise-free  step,  noise-free 
ramp,  noise-free  acceleration,  and  noisy  step  were  applied 
with  the  results  shown  in  Figures  34  through  37  inclusive. 
The  main  point  of  interest  in  these  results  is  that  step 
response  in  the  presence  of  noise  for  this  staleness  factor 
design  is  perhaps  the  best-behaved  so  far,  although  it  is 
not  good  by  any  means.  In  addition,  noisy  inputs  of  ramp 
and  acceleration  were  applied,  but  both  of  these  gave  un¬ 
stable  responses  as  in  previous  designs. 

The  final  design  effort  to  be  described  in  this  work 
was  directed  toward  achieving  ripple-free  response  to  a 
velocity  input  when  sampling  at  a  rate  of  ten  samples  per 
second.  For  this  design  all  of  the  zeros  of  G(z)  must  be 
contained  in  K(z),  so  we  begin  by  writing 
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K(z)  =  (l+5.48614151z~1)(l+.47316625z“1) (1+.05265&21Z”1) 
(l+.ooisggosz"1)^^"1  +  a2z"2) 

Since  we  are  designing  for  zero  steady  state  error  for  a 

ramp  input,  then  it  must  be  true  that 

l-K(z)  =  (l-z"1)2(l  +  b^-1  +  b2z“2  +  b3z“3  +  b4z“4) 

After  multiplying  out  and  equating  like  coefficients  of  z 
we  obtain  six  simultaneous  equation  in  six  unknowns  which 


are  presented  below  in  matrix  form. 


1.0 

0 

1.0 

0 

0 

0 

al 

2.0 

6.013363 

1.0 

-2.0 

1.0 

0 

0 

a2 

-1.0 

2.918063 

6.013363 

1.0 

-2.0 

1.0 

0 

bl 

0 

.140759 

2.918063 

0 

1.0 

-2.0 

1.0 

b2 

0 

.000191 

.140759 

0 

0 

1.0 

-2.0 

b3 

0 

0 

.000191 

0 

0 

0 

1.0 

b4 

0 

— 

_  m 

This  matrix  equation  was  solved  by  computer  using  a  Gaussian 
elimination  library  routine.  With  the  results  obtained  there¬ 
from,  the  coefficients  of  the  process,  K(z),  were  calculated 

as  follows:  A(l)  =  2  -  =  .319530814 

A( 2 )  =  2b,  -  1  -  b2  =  1.70120542 

A( 3 )  =  2b 2  -  b1  -  b3  =  -.39202846 

A( 4 )  =  2b 3  -  b2  -  b4  =  -.59772472 

A( 5 )  =  2b4  -  b3  =  -.03094094 

A( 6 )  =  -b4  =  -.00004212 

These  process  coefficients  were  supplied  to  the  simulator 
program  along  with  the  coefficients  of  P(z)  and  Q(z)  and 
the  values  of  the  four  plant  zeros  to  be  included.  The 
program  then  calculated  D(z)  for  the  single-rate  sampled 
data  system.  For  a  test  input  of  a  noise-free  ramp  the 
output  has  no  overshoot  or  inter-sample  ripple  and  reaches 
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zero  error  in  0.4  seconds.  See  Figure  38.  The  response 
to  a  noise-free  step  is  shown  in  Figure  39 ,  where  the 
"tuned"  nature  of  the  design  is  again  exhibited.  However, 
after  the  initial  large  overshoot  the  step  response,  too, 
is  essentially  ripple  free  and  steady  state  in  0.4  seconds. 

For  the  noisy  ramp  signal  shown  in  Figure  40  the  sam¬ 
pled  data  system  was  unstable,  while  for  a  step  signal  with 
the  same  amount  of  additive  noise  the  response  can  be  seen 
in  Figure  41.  Note  that  the  time  scale  in  Figure  41  is  five 
times  longer  than  in  the  noise-free  experiments,  and  the  re¬ 
sponse  to  the  noisy  signal  shows  no  signs  of  settling  down 
any  further  at  the  end  of  five  seconds  than  it  did  at  the 
end  of  the  first  second.  Clearly,  noise  of  this  power  cannot 
be  tolerated  in  the  input  signal  to  a  sampled  data  system 
regardless  of  the  digital  controller  design. 
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Figure  28 
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Figure  29 
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Figure  30 
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Figure  31 
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Figure  35 
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Figure  36 
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Figure  37 
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Figure  39 
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Figure  40 
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Figure  41 
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8.  Conclusions. 

The  behaviour  of  a  unity  feedback  control  system  he. a 
been  simulated  using  a  digital  computer  program.  Uc mg  this 
simulator  program  studies  have  been  made  of  digital  control 
methods  for  various  plants,  with  particular  attention  fo¬ 
cused  on  the  control  of  a  Naval  gun  mount.  The  azimuth  func¬ 
tion  of  the  gun  mount  is  described  by  a  fifth  order  a  infer¬ 
ential  equation,  and  it  is  believed  that  the  concic  rons 
pertaining  to  it  are  applicable,  at  least  in  part,  to  other 
shipboard  weapon  systems. 

The  use  of  multi-rate  digital  controllers  wa-  inves¬ 
tigated.  In  every'’  instance,  the  double-rate  controller  pro¬ 
vided  faster  response  by  one-half  a  sampling  period.  How¬ 
ever,  only  in  the  special  case  of  ripple-free  dec-  .  t.go  for  a 
step  input  was  the  double-rate  response  an  unalloyed  im¬ 
provement  cvci  the  single-rate  response.  In  genet 
transient  concroj  and  inter-sample  ripple  was  not  uetice- 
ably  improved  using  a  double-rate  controller  fez  the  gun 
mount.  For  other  plants  of  different  frequency  response, 
advantage  mi  gat  oe  taken  of  the  fact  that  inter-sampie 
ripple  when  using  a  double-rate  controller  occurs  ch  t.  ice 
the  error- samp ling  rate. 

Pulse  transfer  functions  (or  Z-transf orrns }  ;  a  func¬ 

tion  of  the  sampling  period.  It  was  found  that  i or  certain 
sampling  rates  valid  simplifying  assumptions  could  made 
which  facilitated  the  design  procedure.  Care  should  oe 
exercised  that  these  simplifying  assumptions  are  net  made  on 
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the  basis  of  a  sampling  rate  which  is  too  slow  in  relation 
to  the  damped  natural  frequency  of  the  system  concerned. 

In  this  regard  the  approximations  made  to  the  pulse  trans¬ 
fer  function  of  the  gun  mount  were  inappropriate  since  they 
were  based  on  a  sampling  rate  roughly  equal  to  the  damped 
natural  frequency  instead  of  six  to  twelve  times  greater  as 
must  be  the  case  for  adequate  control. 

Extreme  sensitivity  of  sampled  data  systems  to  noise 
in  the  input  signal  was  demonstrated.  No  satisfactory  method 
of  coping  with  input  noise  by  design  of  the  over-all  process 
was  found.  So-called  "variance  reduction"  processes  showed 
little  merit  as  the  result  of  these  simulation  studies.  The 
design  criteria  which  came  closest  to  producing  satisfactory 
response  to  a  noisy  step  input  was  the  prototype  with  stale¬ 
ness  factor  of  0.5. 

Without  compensation,  either  analog  or  digital,  the 
gun  mount  is  closed-loop  unstable.  Several  digital  control 
functions  were  developed  which  gave  stable  response  to  noise- 
free  input  signals.  None  of  these  gave  stable  operation  when 
noise  was  present  in  a  first  order  (ramp)  input  or  a  second 
order  (acceleration)  input,  but  for  a  noisy  zero  order  (step) 
input  digital  control  achieved  a  response  that  was  at  least 
bounded  for  the  duration  of  the  experiment,  albeit  unsatis¬ 
factory  for  most  conceivable  applications.  It  is  concluded 
that  sampled  data  system  sensitivity  to  additive  noise  in  the 
input  increases  with  the  order  of  the  input  signal. 

The  computer  program  described  in  Appendix  II  for 
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finding  Z-transf orms  was  very  helpful  in  this  investigation. 
It  is  hoped  that  this  program,  as  well  as  the  control  system 
simulator  upon  the  results  of  which  this  work  was  based, 
will  prove  useful  to  others  in  the  field  of  sampled  data 
systems.  It  is  further  hoped  that  this  work  has,  in  a  small 
way,  indicated  the  truthfulness  of  the  assertion  that  digital 
compensation  techniques  introduce  a  degree  of  flexibility  and 
accuracy  in  a  control  system  that  is  not  presently  achievable 
with  state-of-the-art  analog  techniques. 
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APPENDIX  I 


DIGITAL  COMPUTER  PROGRAM  FOR  SIMULATION  OF  A  SAMPLED 

DATA  SYSTEM 


General  Description  of  Program  ALSTAPP 

This  program  is  a  simulation  of  the  hybrid  control 
system  shown  below. 


-3 


K(S) 


Three  choices  are  available  for  the  test  input?  namely , 
a  unit  step,  a  unit  ramp,  and  a  unit  acceleration.  The 
plant  may  be  up  to  thirtieth  order  and  is  currently  sub¬ 
ject  to  the  following  restrictions: 

(1)  No  plant  zeros  are  allowed,  i.e.,  the  numerator 
of  H(s)  must  be  a  constant. 

(2)  The  plant  may  not  have  a  pole  at  the  origin  higher 
than  first  order,  i.e.,  limited  to  Type  0  and  Type  1  servos. 

(3)  Only  the  constant  term  may  be  missing  in  the  dif¬ 
ferential  equation  describing  the  plant. 

In  essence,  the  plant  is  restricted  to  having  a  trans¬ 
fer  function  of  the  general  form 


K(s)  = 


VK 


s(s+  pkxs+et)~‘(s+  Pk-() 


simple  poles.  Pi 


H  (s)  =  A(K+l)  SK  +  A(k)  sk_1  +  •  •  •  t  >4(2)  s  +■  A  (l) 

where  A(l)  is  normally  zero  for  Type  1  servo  and  is  the  only 

missing  term  allowed. 
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The  rate  at  which  the  raw  error,  El,  is  sampled  can 
be  selected  at  will0  This  sample  period  T  is  chosen  by 
the  user  and  supplied  to  the  computer. 

DELT  is  the  small  increment  of  time  used  in  the  simu- 
lation  of  the  continuous  plant  by  numerical  solution  of  the 
plant's  differential  equation.  This  is  done  one  thousand 
times,  so  that  the  total  problem  time  is  1000  x  DELT. 

DELT  should  be  chosen  as  small  as  possible  consistent  with 
the  desired  total  problem  time. 

Although  one  thousand  values  of  the  solution  are  gen~ 
erated  and  are  available  for  use  by  a  computer  plotting 
routine,  only  every  tenth  solution  is  produced  on  the  Ana- 
lex  printer,  since  it  was  felt  that  a  column  of  one  hundred 
values  is  sufficient  for  perusal  by  eye  and  more  than  ade¬ 
quate  for  plotting  by  hand. 

This  program  performs  four  major  tasks.  The  first 
three  are  performed  every  time  the  program  is  executed, 
but  it  is  optional  with  the  user  whether  or  not  the  fourth 
task  will  be  done.  These  tasks  are: 

(1)  Simulation  of  a  continuous,  un-compensated,  unity 
feedback  system  for  the  given  plant. 

(2)  Calculation  of  coefficients  for  a  single-rate 
digital  controller. 

(3)  Simulation  of  the  sampled-data  system  employing 
a  single-rate  digital  controller  using  the  coefficients 
calculated  above. 

(4)  Simulation,  if  desired,  of  a  sampled-data  system 
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employing  a  double-rate  controller.  The  coefficients  of 
the  double-rate  controller  must  be  supplied  by  the  user,, 
and  if  they  are  not,  the  fourth  task  is  not  performed. 

Information  must  be  supplied  to  the  program  via  data 
cards.  Detailed  instructions  for  the  preparation  and  con¬ 
tents  of  these  cards  is  given  later.  Seventeen  (17)  cards 
are  required  if  the  program  is  to  do  all  four  tasks*  while 
only  the  first  thirteen  (13)  need  be  supplied  if  simula¬ 
tion  of  the  double  rate  controller  is  not  desired.  CAUTION  § 
The  cards  must  be  presented  to  the  computer  in  exactly  the 
order  and  form  specified.  The  slightest  error  will  cause 
faulty  operation  or,  perhaps,  no  operation  at  all. 

Discussion  of  Information  Supplied  by  the  User 

The  general  nature  of  the  information  that  must  be 
supplied  to  the  computer  is  described  next. 

(1)  The  user  must  choose  the  small  time  increment 
DELT  to  be  used  in  the  numerical  method  solution  of  the 
plant's  differential  equation.  The  smaller  this  increment 
the  more  accurate  the  solution.  However,  the  total  time 
simulated  will  be  only  one  thousand  times  DELT,  s©  the  user 
should  first  determine  the  time  span  of  interest. . „ 1. e. , 

an  estimate  ©f  how  long  it  will  take  the  system  to  achieve 
steady  state,  or  how  much  of  the  transient  condition  it  is 
desired  to  observe  for  lightly  damped  systems. . „ then  divide 
this  total  time  by  1000  to  obtain  the  value  for  DELT. 

(2)  The  user  must  choose  the  sampling  period  to  be 
employed  in  the  sampled-data  system  with  single-rate 
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controller  simulation.  This  value  is  called  T. 

(3)  The  user  chooses  the  test  input  to  be  applied 

to  the  system.  Unit  values  of  step,  ramp,  and  acceleration 
are  available. 

{NOTE  Is  These  first  three  pieces  of  information,  namely, 
DELT,  T,  and  INPUT  are  supplied  by  data  card  number  one. ) 

(4)  The  specifications  of  the  continuous  plant  are 
supplied  to  the  computer  by  data  cards  2,  3,  and  4.  The 
order  of  the  plant  is  denoted  by  K.  This  is  the  highest 
power  of  the  Laplace  variable,  S,  in  the  denominator  of 
the  plant  transfer  function,  H(s).  The  numerator  gain 
constant  of  H(s)  is  denoted  by  VK.  The  order  of  the  pole 
at  the  origin  is  denoted  by  NTYPE  and  is  limited  to  0  or  1. 
These  values,  K,  VK,  and  NTYPE  are  supplied  by  data  card 
number  two. 

(5)  H(s)  must  be  known  in  the  polynomial  form,  that 

is  to  say,  a  numerator  gain  constant,  VK,  over  a  denominator 
polynomial  in  S,  The  constant  coefficients  of  this  denomi¬ 
nator  polynomial  are  denoted  by  A(i)  and  must  be  supplied 
to  the  computer.  The  only  coefficient  which  may  have  a 
zero  value  is  the  constant  term  A(l).  The  number  of  co¬ 
efficients  is  denoted  by  NA  and  is  supplied  by  data  card 
number  3.  Even  though  A(l)  is  normally  zero  (in  a  Type  1 
servo)  it  must  be  included  in  the  count,  NA,  on  card  3  and 
its  zero  value  shown  in  the  first  ten  columns  of  data  card 
number  four  (4).  In  the  usual  case  of  a  iype  1  servo,  the 
value  of  NA  is  one  more  than  the  order,  K,  of  the  plant. 
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(6)  For  the  simulation  of  the  sampled-data  systems,, 

Z-transform  techniques  are  implemented*  Data  reconstruc- 

tion  is  accomplished  by  a  simulated  zero  order  hold*  which 

must  be  included  in  the  Z-transform  approach  to  the  plant* 

If  we  call  the  plant  and  hold-circuit  combination  G(s),  then 
-  sT 

\  “  €  and  the  pulse  transfer  func- 

s 

tion  (or  z-transform)  is 


«...  Z  T  ■  8M 

This  pulse  transfer  function  for  the  plant  and  hold  must 

be  supplied  by  the  user.  Its  general  form  is  a  ratio  of 

-1 

polynomials  in  Z  as  shown  below: 

“  MP 

G(z)  =  ■>■•••  P(Nf)  z _ 

’  0(a)  ~  Q(t)  +  +Q(3)£-2-+...^ 

where  NP  is  the  number  of  terms  (and  also  the  highest 
power  of  7T  )  in  the  numerator,  P(Z),  and  NQ  is  the  num¬ 
ber  of  terms  in  the  denominator,  Q(Z).  In  order  for  G(Z) 
to  be  a  physically  realizable  pulse  transfer  function,  the 
term  Q(l)  must  be  a  constant,  and  in  fact,  is  usually  unity. 

(7)  Although  G(Z)  is  supplied  to  the  computer  in 
polynomial  form,  it  must  be  known  in  factored  form,  in 
order  that  any  poles  or  zeros  of  G(Z)  that  lie  on  or  out¬ 
side  the  unit  circle  in  the  Z-plane  may  be  accounted  for. 
Usually,  the  plant  is  open  loop  stable  and  only  zeros  need 
be  accounted  for,  but  the  program  will  also  process  infor¬ 
mation  and  provide  the  digital  controller  coefficients 
when  the  plant  has  a  pole  outside  the  unit  circle.  All 
such  zeros  must  be  supplied  by  the  user  under  the  generic 
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term  ZERO(i)  and  the  number  of  them  is  denoted  by  NZ. 

Up  to  ninety  “-nine  poles  and/or  zeros  may  be  accounted  for* 
(8)  The  overall  pulse  transfer  function  of  the  en¬ 
tire  system  is  denoted  by  K(Z)  and  is  defined  as 

CO) 


K(z)  * 


R(a) 


For  the  purposes  of  this  program,  K(Z)  must  be  a 
numerator  polynomial  only.  This  implies  a  "finite-memory" 
process.  If  an  infinite-memory  or  "non-finite  settling 
time  process"  is  desired,  then  it  must  be  approximated 
by  a  finite  number  of  terms  of  the  infinite  series  ob¬ 
tained  by  long  division.  The  program  will  operate  on  up 
to  ninety-nine  (99)  such  terms  for  K(Z).  The  general  form 
for  K(Z)  is 

—  I  "  N  K 

«(2)  -  AK(i)  2  +  AK(a)  2  ■+•***+  AK  (nk)  2 

where  NK  is  the  number  of  terms  provided.  K(Z)  must  be 
supplied  by  the  user. 

NOTE  2:  In  order  to  supply  an  array  of  numerical  data  to 
the  computer  two  (2)  data  cards  are  required.  The  first 
card  tells  the  computer  how  many  values  are  in  the  array, 
and  the  next  card  (or  cards)  supplies  these  values.  As¬ 
suming  for  the  moment  that  no  array  contains  more  than 
eight  (8)  values,  the  following  information  is  supplied 
to  the  computer  by  the  data  cards  listed  below s 

Card  Nos.  Information 

3  &  4  Array  of  coefficients  of  plant 

differential  equation  (i.e. , 
denominator  of  H( s ) ) . 


5  &  6 


Array  of  coefficients  of  nu¬ 
merator  of  G(Z)0 
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Information 


Card  Nos, 
7  &  8 


9  &  10 

11  &  12 


Array  of  values  of  zeros  and 
poles  of  G(Z)  to  be  included 
in  overall  process  pulse  trans- 
fer  function,  K(z), 

Array  of  coefficients  of  de= 
nominator  of  G(Z)e 

Array  of  coefficients  of  over= 
all  process  pulse  transfer 
function,  K(Z), 


(9)  For  certain  investigations  it  may  be  desirable 
to  limit  the  output  of  the  digital  controller  t©  a  signal 
of  a  given  magnitude  to  avoid  plant  saturation,,  If  so, 
the  user  must  supply  this  value  as  a  positive  decimal 
value,  E2KAX,  and  the  program  will  limit  the  processed 
error  from  the  controller  to  plus  or  minus  this  value. 

If  no  limiting  is  desired,  a  value  of  zero  (0*0)  should 
be  given  to  E2I4AX  and  the  program  then  sets  a  magnitude 
limit  of  one  million  on  the  processed  error.  Data  Card 
13  contains  the  data  on  E2MAX. 

NOTE  3:  In  very  simple  cases  it  may  be  that  there  are  no 
zeros  or  poles  of  G(Z)  that  need  to  be  included  in  K(Z)0 
When  this  occurs,  obviously  the  datum  NZ  is  zero  and  the 
array  ZERO(i)  does  not  exist.  Under  these  conditions, 
prepare  data  card  7  as  shown  below 


C?,l,  1. 


Card  7 


00 


and  remove  card  8  completely! 
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Card  No 


DATA  CARDS  for  Input  to  Program  ALSTAPP 


Format  Contents  Remarks 


1 

2F10o0#Il 

DELT,  T,  INPUT 

Total  time  ■  1000 

DELT0/  Sample  period 

T  ss  ly/ error  sampling 
rate©  DELT  is  the  in- 
crement  used  for  simu¬ 
lation  of  continuous 
plant.  INPUT  *  1  for 
a  step,  2  for  a  ramp, 

3  for  acceleration 0 

2 

I2,F18.0,I1 

K,  VK,  NTYPE 

K  =  order  of  pi ant e 

VK  -  plant  gain© 

NTYPE  s  type  ©£  servo 

C i © © o ,  1 ,  ^  @  ^ ) 

3 

12 

NA 

NA  “  no©  terms  in  de- 
nominator  of  plant 
transfer  function 8 
Usually  NA  -  K  +  i0 

H(S)  a 

VK 

A(K+1)SK  +0. „A(2)S+A(1) 

4 

8F10C  0 

A( I )  I  =  1,  NA 

A(I)  is  the  array  of 
coefficients  ©f  plant 
differential  equation 
denominator  of 
H(S)  above)© 

5 

12 

NP 

NP  -  no©  terms  in  nu¬ 
merator  ©f  Z-transform 
of  plant  and  zero- 
order  hold  whose  z= 
transform  is 

G(z\  -  P(.Z Jt  _ 

G{z}  “  off! 

6 

8FX0.0 

P(I)  1=1,  NP 

PU)z4+P{2)z“2+oe. 

1+Q(2)z“1"JQ(3)z“2+bo. 

+  p(np)z~np 

+Q(NQ) z-NQ+l 

7 

12 

NZ 

NZ  =  noe  zeros  of  G(z) 
to  be  included  in  over¬ 
all  pulse  transfer  fen, 
K(Z). 
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DATA  CARDS  for  input  to  Program  ALSTAPP  ( Continued) 


Card  No. 

Format 

Contents 

Remarks 

8 

8F10.0 

ZERO(I)  1=1, NZ 

ZERO(I)  =  array  of 
values  of  such  zeros 
of  G(Z)  as  are  to  be 
included  in  K(Z). 

9 

12 

NQ 

NQ  =  no.  terms  in  de¬ 
nominator  of  G(Z). 

10 

8F10.0 

Q ( I )  1=1,  NQ 

Q( I )  =  array  of  de¬ 
nominator  coefficients. 

11 

12 

NK 

NK  =  no.  terms  in  over- 
all  pulse  transfer 
function 

K(Z)=AK(l)z“1+AK(2)z“2 
.. .+AK(NK)z“NK 

12 

8F10.0 

AK(  I )  I  =  1,  NK 

AK( I )  -  array  of  pro¬ 
cess  coefficients 
shown  above. 

13 

F10.0 

E2MAX 

E2MAX  is  the  limit 
placed  on  the  absolute 
value  of  the  output  of 
the  digital  controller. 

Notes  to  the  User: 

The  FORMAT  conveys  information  concerning  the  mode  and 
location  on  the  IBM  card  for  each  item  of  numerical  infor¬ 
mation  to  be  supplied  to  the  program.  If  the  following 
remarks  are  not  sufficient,  refer  to  the  CDC  Fortran  Manual. 

"I"  means  an  integer  value.  Whole  numbers  only.  Do 
not  show  a  decimal  point.  The  integer  value  must  be  "right 
justified"  in  its  allotted  field.  If  the  field  is  left 
blank,  a  value  of  zero  will  be  inferred. 

"F"  means  a  decimal  value,  either  less  than  or  greater 
than  unity.  A  decimal  point  must  be  shown  somewhere  in  the 


103 


allotted  field*,  If  the  value  is  negative,  the  minus  sign 
must  also  fall  within  the  field.  Precise  location  of  the 
value  (including  sign,  if  any,  and  the  decimal  point  al¬ 
ways)  within  the  field  does  not  matter. 

The  numbers  following  I  and  F  give  the  number  of  col- 
umns  on  an  IBM  card  that  are  reserved  for  the  field  of 
that  particular  datum.  For  example,  take  the  FORMAT  con¬ 
trol  for  data  card  number  2:  Columns  1  thru  2  are  reserved 
for  the  integer  value  of  K.  If  the  integer  value  of  K  is 
four,  then  card  2  should  have  04  in  its  first  two  columns. 
The  next  eighteen  columns  (3  thru  20)  on  card  number  2 
must  contain  1500000.0  somewhere  between  columns  3  and  20 
inclusive.  On  both  cards  1  and  2,  column  21  must  contain 
a  single  integer. 

On  those  cards  having  an  12  format,  columns  1  and  2 

contain  an  integer  value  between  one  and  ninety-nine.  If 
the  value  is  less  than  ten,  it  must  appear  in  column  2 
(i.e.,  the  integer  value  must  be  right  justified  in  its 
field.  Decimal  values  need  not  be  justified).  In  every 
case  these  12  integer  values  tell  the  computer  how  many 
decimal  values  t©  read  from  the  next  card  or  cards.  This 
brings  us  t©  the  FORMAT  8F10.0.  As  before,  the  F1Q.0 
means  a  decimal  value  with  ten  card  columns  reserved  for 
it.  8F10.0  means  there  may  be  up  to  eight  such  values  on 
a  single  card  (a  card  has  eighty  columns),  and  the  com¬ 
puter  will  read  another  card  if  it  has  to,  in  order  to 
read  as  many  values  as  were  specified  by  the  previous  12 
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card.  Most  input  arrays  will  contain  eight  or  less  values, 
however,  and  thus  will  fit  nicely  on  a  single  card.  Note 
that  these  decimal  fields  of  ten  columns  each  begin  in 
card  columns  1,  11,  21,  31,  ...etc.  Card  1  has  two  dec±~ 
mal  fields  of  ten  columns  each,  followed  by  an  integer 
field  of  one  (1)  in  column  21.  Card  13  has  a  single 
decimal  field  occupying  card  columns  1  thru  10. 

If  a  decimal  value  of  zero  is  to  be  entered,  then 
place  0.0  in  the  appropriate  field.  For  example,  suppose 
array  A  has  three  decimal  values  A(l)  -  zero,  A{2)  * 
twenty,  A(3)  -  minus  two  and  a  half.  Then  Card  4  should 
contain 

Between  Columns  The  Decimal  Value 

1  and  10  0.0 

11  and  20  20.0 

21  and  30  -2.5 

If  no  sign  is  shown,  the  value  is  taken  as  being 

positive. 

In  order  to  use  this  program  to  test  a  system  using 
a  single- rate  controller,  only  the  first  thirteen  data 
cards  are  required. 

If  you  desire  to  test  a  system  using  a  double-rate 
controller,  then  an  additional  four  data  cards  are  re¬ 
quired.  These  cards  contain  information  as  to  the  double¬ 
rate  controller  constants  and  follow  card  13  (no  inter¬ 
vening  blanks  8 )  as  shown  below? 
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ADDITIONAL  DATA  CARDS  FOR  DOUBLE-RATE  CONTROLLER 


Card  No.  Format  Contents  Sana rks 


14 

12 

NAA 

NAA  is  the  no.  terms 
in  numerator  of  double¬ 
rate  controller,  D(Z^). 

15 

8F10.0 

AA( I )  1=1 ,  NAA 

Array  of  numerator 
coefficients. 

16 

12 

NBB 

NBB  is  the  no.  terms 
in  denominator  of 

D(z 

2  BB(X)*BB{2)z“2+ 

.. .+aa(naa)z“naa+1 
...SB(Na3)z-NBE+l 

17 

8F10.0 

BB( I )  1=1,  NBB 

Array  of  denominator 
coefficients. 

Example  Problem  usinq  Program  ALSTAPPg 

l  |  i  item  ,  ft  - - —  ■  I—  ■  i  fft  I  I  -i  ffi»  i  J  ■  i  ■i»A— —  minim  ■■hi  ir  — i 


H(S) 


80, 000 

S3  +  95S2  +  1750S 


3y  looking  at  H(s)  we  can  see  that  VK  =  80G00.0 

K  **  3 
NA  *  4 
A(l)  »  0.0 

A(  2 )  *  1750.0 
A(  3 )  *  95.0 
A(4)  a  1.0 
NTYPE  -  1 


Let  us  choose  to  observe  the  response  of  this  plant  for 
one  second  after  the  application  of  a  unit  step  Input. 
So  DELT  ss  1/1000  -  0.001 
INPUT  »  1 
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Let  us  select  a  sampling  rate  of  100  times  per  sec¬ 


ond  for  the  sampled  data  system,  i.e. ,  T  *  QoCl., 

We  prepare  the  first  four  data  cards  as  shown  below* 


Column 

1 


Column 

11 


Column 

21 


31 


By  partial  fraction  expansion  and  tables  of  Z-tr ans- 
forms  we  find  the  pulse  transfer  function  of  the  plant- 

and-hold  to  ha 

in  factored  forms 

i  i  + 


G(  Z )  = 

(1“S  : 


209327Z  )  (1  +  20  96727 3Z 

^0i778801Z~1)  (1 


0o 49658 5Z  l) 


in  polynom Lai  foritu 

G(7)  =  «  010600Z°;1  +  „  033672Z~2  +  .006584Z"3  _  . 

1.0  -  2.2753862“*  +  1.662127Z*2  -  0o 3867412“ 3 

There  is  one  aero  outside  the  unit  circle  at  -2*967275. 

Choosing  a  minimal  prototype  process  we  find  that 

K(Z)  “  0*2520S2Z~1  +  0o747938Z“2 

and  we  see  the  following  data: 

NP  -  03  P(l)  =  .010600 

P( 2 )  =  .033672 
P( 3 )  =  .006584 

NZ  s  0.1  ZERO  (1)  =  -2.967275 
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NQ  =  04 


NK  -  02 

So  we  prepare 


Q(l) 
Q(  2) 
Q(  3) 
Q(4) 

AK(1) 

AK(2) 

data  cards  as 


=  1.0 

=  -2.275386 
=  1.662127 
=  -0.386741 

=  0.252062 
=  0.747938 

follows? 


Column  Column  Column  Column 

1  11  21  31 


Card  12 
Card  11 
Card  10 
Card  9 
Card  8 
Card  7 
Card  6 
Card  5 


Suppose  we  do  not  desire  to  limit  the  output  of  the 

digital  controller.  Then  set  E2MAX  =  zero  by  preparing 

card  13  as  shown  below? 

Column 

1 

/  (0.0  Card  13 

( 

If  it  is  not  desired  to  simulate  a  double-rate  digital 
controller *  then  only  the  first  thirteen  (13)  data  cards 
are  required. 

If  double-rate  action  is  desired*  then  the  user  must 
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supply  the  coefficients  of  the  double-rate  controller. 

The  pulse  transfer  function  of  the  double-rate  controller 
has  the  following  general  form: 

D(Z0)  =  AA(1*  +  AMSJZ"1  +  AA(3)Z“2  +  ...  +  AA(NAA)  ZjNAA+1 
BB(  1 )  +  BB(2)Zj1  +  BBOJZ^2  +  ...  +  B3(NBB)Z“NBB+1 

where  all  coefficients  have  been  normalized  with  respect 
to  BB( 1 )  and  BB(1)  -  1.0. 

The  format  for  cards  14  thru  17  is  identical  to  that 
for  cards  5  thru  12. 

In  the  double-rate  controller  the  error-sampling 
period,  T,  remains  the  same  as  it  was  for  the  single-rate 
case.  It  is  the  output  from  the  controller  that  occurs 
at  the  double-rate. 

Modification  to  allow  for  addition  of  white  noise  to  the 
test  input. 

When  used  as  described  in  the  main  body  of  this 
appendix,  the  program  simulates  a  control  system  whose 
input  is  completely  deterministic.  Gaussian  noise  of 
selected  mean  and  variance  may  be  added  to  the  input  by 
making  an  addition  to  the  first  data  card  and  inserting 
a  new  data  card  between  the  previous  card  1  and  20  as 
follows. 

1.  The  presence  of  noise  in  the  input  must  be  re¬ 
quested  by  placing  the  integer  1  in  column  31  of  the  first 
data  card  (which  already  contains  DELT,  T,  and  INPUT). 

2.  The  new  data  card,  which  must  follow  card  1 

above,  contains  an  integer  number  between  01  and  99  right 
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justified  in  columns  1  and  2,  and  a  decimal  value  in 
columns  3  through  20.  The  decimal  value  specifies  the 
selected  mean  value  of  the  noise,  and  in  most  cases  will 
be  shown  as  0.0. 

The  integer  number  specifies  the  number  of  uniform- 
ly  distributed  random  numbers  from  which  each  normally 
distributed  random  number  is  obtained.  For  a  good  approx- 
imation  to  white  noise,  this  number  should  be  between  ten 
and  twenty.  The  variance  of  the  Gaussian  distribution 
depends  on  this  number,  N,  as  shown  below? 

Variance  -  ^ 

As  N  approaches  one,  the  noise  approaches  a  uniform 
distribution  between  -1/2  and  +1/2,  for  a  selected  mean  of 

zero. 

The  remainder  of  the  data  cards  are  then  included 
in  exactly  the  order  and  format  described  earlier. 

Modification  to  allow  the  program  to  use  single-rate  digi¬ 
tal  controller  coefficients  supplied  by  the  user. 

When  used  as  described  in  the  main  body  ©f  this 
appendix,  the  program  computes  the  coefficients  of  the 
single-rate  controller  using  the  following  formula? 

0£Z)  -  .  !  2  kTz)' 

In  order  to  do  this,  the  user  must  supply  data  on 
the  numerator  and  denominator  of  G(Z),  the  overall  pulse 
transfer  function,  K{Z),  and  which  zeros,  if  any,  of  G(Z) 
that  are  included  in  K(Z)„  CAUTION:  Slide-rule  accuracy 
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is  usually  not  sufficient  except  in  very  simple  cases i 
It  is  recommended  that  a  desk  calculator  be  used  to  pro¬ 
vide  this  data  accurately  to  at  least  five  decimal  places. 

In  case  the  user  wishes  to  calculate  his  own  con¬ 
troller  coefficients  for  the  single-rate  sampled-data 
system  simulation^  he  may  cause  the  program  to  omit  the 
calculation  of  D(Z)  and  to  employ  the  controller  supplied 
by  the  user.  In  order  to  do  this,,  the  data  cards  are  ar¬ 
ranged  as  shown  below? 


rd  No. 

Contents 

Remarks 

1 

DELT,  T,  INPUT 

Same  as  before 

2 

K,  VK 0  NTYPE 

Same  as  before 

3 

NA 

Same  as  before 

4 

A(I)  I  =  1,  NA 

Same  as  before 

5 

Blank  card! 

Indicates  D(Z)  supplied  by 
user 

6 

NAA 

Single-rate  controller 
efficients 

CO- 

7 

AA( I )  1=1,  NAA 

Single-rate  controller 
efficients 

CO- 

8 

NBB 

Single-rate  controller 
efficients 

co¬ 

9 

BB( I )  I  -  1,  NBB 

Single-rate  controller 
efficients 

co¬ 

10 

32MAX 

Same  as  before 

11* 

NAA 

Double-rate  controller 
efficients 

co- 

12* 

AA(  I )  I  ■-  1,  NAA 

Double-rate  controller 
efficients 

CO- 

13* 

NBB 

Double-rate  controller 
efficients 

CO- 

Ill 


Card  No. 


Contents 


Remarks 


14*  3B(I)  I  =  1,  N3B  Double-rate  controller  co- 

off  icients 

♦Optional. 

The  first  four  cards  are  just  the  same  as  they  would 
be  for  a  deterministic  input  where  D(Z)  is  to  be  computed 
by  the  program.  A  blank  fifth  card  signals  the  program 
that  Z>(Z)  will  be  supplied  by  the  next  four  cards.  These 
next  four  cards  (6,  1,  Q,  and  9)  give  the  data  for  tne 
single-rate  controller  exactly  as  described  earlier  in  the 
main  body  of  this  memo  for  the  double-rate  controller  co¬ 
efficients.  Card  10  gives  3 2 MAX.  If  double-rate  action 
is  also  desired,  the  last  four  cards  must  give  the  double¬ 
rate  D(Z). 
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FLOW  DIAGRAM  for  Program  ALSTAPP 
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APPENDIX  II 


DIGITAL  COMPUTER  PROGRAM  FOR  FINDING  Z-TRANSFORMS  AND 

ROOTS  OF  POLYNOMIALS 

General  Description  of  Program  Stappl 

This  program  is  a  combination  rootfinder  and  Z-trans- 
form  taker.  It  will  operate  in  one  of  three  modes ,  as 
specified  by  the  user.  These  modes  are: 

Mode  1 

Program  operates  as  a  rootfinder  only.  The  user  sup¬ 
plies  the  degree  and  coefficients  of  the  polynomial  whose 
roots  are  to  be  found.  Handles  polynomials  up  to  and  in¬ 
cluding  thirtieth  degree. 

Mode  2 

It  is  desired  to  take  the  Z-transform  of  a  transfer 

function  which  is  not  known  in  factored  form,  such  as 

H(s)  =  _ VK _ _ 

A(N+1)sN+A(N)sN"1+  .  .  .  +  A( 2) s  +  A( 1 ) 

The  program  first  finds  the  roots  of  the  denominator 

polynomial.  Then  it  takes  the  Z-transform  of  the  above 

plant  transfer  function  preceded  by  a  zero-order  hold,  i.e. 

G<z>  =  &  (¥'st-H(s))  -  =  §{§} 
and  presents  the  z-plane  poles  of  G(Z)  as  well  as  the 
numerator  and  denominator  polynomials  in  Z”^.  Then  it 
operates  as  a  rootfinder  on  the  numerator,  P(Z),  in  order 
to  display  the  zeros  of  G(Z). 

Mode  3 

It  is  desired  to  find  the  Z-transform  of  a  plant  and 
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hold  combination  when  the  transfer  function  of  the  plant 
is  known  in  factored  form,  such  as 


As  in  Mode  2,  the  Z-transform  is  taken  with  the  plant 
preceded  by  a  zero-order  hold,  i.e. 


The  program  displays  the  Z-plane  poles  of  G(Z)  and  the 
numerator  and  denominator  polynomials.  It  then  operates  as 
a  rootfinder  on  the  numerator  polynomial,  P(Z),  in  order  to 
display  the  zeros  of  G(Z). 

Restrictions. 

Mode  1  is  relatively  unrestricted.  It  will  find  the 
roots  (real  and  imaginary  parts)  of  polynomials  in  the 
following  form  (for  N  up  to  and  including  30): 

P(x)  =  A(N+1)  xN  +  A(H]xN=1  +  ...  +A  ( 2 )  x  +  A(  1 ) 

Modes  2  and  3  are  restricted  to  polynomials  of  degree 
ten  or  less,  with  all  roots  lying  on  the  real  axis  and  no 
repeated  roots.  Mode  2  will  accept  a  polynomial  whose 
eleventh  root  is  at  the  origin.  In  both  modes  the  Z- 
transform  operation  is  predicated  upon  the  plant  being  a 
type  1  servo  (i.e.  the  Laplace  transfer  function  has  a 
single  pole  at  the  origin). 

Since  the  program  will  not  take  the  Z-transforrn 
when  roots  are  complex,  Mode  2  checks  the  magnitude  of 
each  imaginary  part  with  respect  to  the  corresponding  real 
part  of  each  root  it  finds.  If  the  magnitude  of  any 
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imaginary  part  exceeds  l/100th  of  the  corresponding  real 
part,  the  program  stops.  If  all  imaginary  parts  of  the 
roots  found  in  the  first  phase  of  Mode  2  are  less  than 
l/100th  the  magnitude  of  their  corresponding  real  parts, 
then  they  are  considered  negligible  and  the  root  is  treat¬ 
ed  as  entirely  real.  The  program  then  proceeds  as  de¬ 
scribed  earlier. 

If#  thru  error,  Mode  2  is  selected  when  the  unfactored 
denominator  of  H(s)  is  higher  than  degree  ten  (exclusive 
of  root  at  origin)  then  the  program  stops  after  finding 
the  roots  of  the  denominator. 

Mathematical  Methods:  Root-finder 

The  rootfinding  techniques  is  a  modification  of  Bair- 

stow's  iterative  procedure  for  finding  a  quadratic  divisor 
2 

of  the  form  x  -PP  x-QQ.  When  this  divisor  is  found  suf¬ 
ficiently  closely,  it  is  solved  by  the  quadratic  formula, 
and  divided  out  to  leave  a  reduced  polynomial  two  degrees 
lower.  The  criterion  of  suitable  closeness  is  applied 
first  to  the  change,  DPP,  in  PP  and  then,  if  PP  is  suf¬ 
ficiently  accurate,  to  the  change,  DQQ,  in  QQD  In  both 
cases  it  is  considered  sufficiently  accurate  when  the 
magnitude  of  the  next  correction  to  the  coefficient  is 
less  than  10”^  of  the  magnitude  of  the  coefficient  itself. 

The  essence  of  Bairstow's  iterative  procedure  is  to 
divide  a  polynomial,  say  A(x),  by  the  trial  quadratic 
divisor  to  obtain  a  first  quotient  polynomial,  say  Q(x), 
plus  a  first  remainder.  Then  this  first  quotient  is 
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divided  again  by  the  same  trial  divisor  to  obtain  a  sec¬ 
ond  quotient,  say  T(x),  plus  a  second  remainder.  Then 
the  corrections  necessary  in  the  coefficients  of  the 
trial  quadratic  divisor  are  determined  from  a  suitable 
combination  of  the  coefficients  of  the  first  and  second 
remainders.  A  symbolic  example  will  illustrate  the  method. 

Let  F(x)  =  A(N+1)xN+A(N)xN“‘*‘+  . . . +A( 2 )x+A( 1 )  =  the  poly¬ 
nomial  whose  roots  are  to  be  found. 

2 

Let  the  trial  quadratic  divisor  be  TD(x)  =  x  -PP.x-QQ 

NOW  F(x)  =  TD(x)« rQ(N+lxN=1+Q(N)xN"2+c . 0+Q(4)x+Q(3)l 
+Q( 2 )x+Q( 1 )  L 

where  the  first  quotient  is  Q(x)  inside  the  brackets  and 
the  first  remainder  is  Q(2)x+Q(l). 

Dividing  again  by  the  trial  quadratic  divisors 

Q(x)  =  TD(x)«  [t(N-1)xN_3+T(N-2)xN"4+...+T(4)x+T(3)1 
+T( 2 )x+T( 1 ) 

where  the  second  remainder  is  seen  to  be  T(2)x+T(l) 

The  corrections  to  the  coefficients  of  TD(Z)  are  computed 
from  the  following  formulae: 

M  =  PP.T(l)  +  QQ. T( 2 ) 

D  =  T(1 )  2-M.  T(  2 ) 

Dpp  =  T( 2) «Q(1)  -  T(  1 ) *Q( 2) 

DQq  =  MfQ(2)  -  T(  1 )  »Q.(1|) 

The  modification  of  Bairstow°s  procedure  came  about  as 
the  result  of  several  runs  using  test  polynomials.  It 
was  found  that  convergence  in  the  case  of  non-repeated 
roots  was  more  consistently  obtained  if  each  set  of  iter¬ 
ations  on  successive  reduced  polynomials  was  begun  from 
the  origin.  In  case  there  are  repeated  roots  the  return 
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to  the  origin  is  delayed  until  all  roots  at  that  location 
have  been  found.  It  was  also  found  that  an  occasional 
set  of  iterations  would  not  terminate  even  though  no  fur- 
ther  refinements  were  being  made  in  the  coefficients  of 
the  trial  divisor.  When  this  occurs  the  program  assumes 
the  trial  divisor  is  satisfactory  after  fifty  iterations. 
No  degradation  of  accuracy  has  been  observed  due  to  this 
artificial  convergence  device. 

Because  of  the  sharpness  of  the  criterion  for  a 
satisfactory  trial  divisor  it  appears  that  accuracy  of 
root  values  to  at  least  seven  significant  figures  can  be 
expected.  Since  the  sum  of  the  roots  of  an  Nth  degree 
polynomial  is  the  coefficient  of  the  term  containing 
the  variable  to  the  (N-l)th  power,,  the  program  calcu- 
lates  a  check  Siam  from  the  root  values  found  and  dis¬ 
plays  this  check  sum  along  with  the  true  value  as  an  in¬ 
dication  of  the  accuracy  achieved.  This  also  provides 
the  user  with  the  means  to  tell  at  a  glance  whether  or 
not  the  procedure  converged. 


Mathematical  Methods:  Z-transformer 

The  Z-transform  portion  of  the  program  is  restricted 
to  type  1  servos  preceded  by  a  zero  order  hold.  That  is, 
the  general  foarm  of  the  function  whose  Z-transform  is  to 
be  taken  must  be 

-sT 

n,  x  _  1-e  S'  .  _ VK 

“  S  s(  s-Px )  (  s-P2 ) . . .  (  s-PN ) 

Where  s  is  the  Laplace  variable ,  T  is  the  sampling  period. 
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and  N  is  the  number  of  poles  on  the  real  axis?  not  in¬ 
cluding  the  pole  at  the  origin. 

Note: 

All  poles  must  be  real  and  none  can  be  repeated. 

sT 

By  definition,  Z=e  ,  so  by  partial  fraction  expansion  we 
can  write 

G(Z)  =  VK>(Z-l)y/  A  B  Cyj.  C£2l  .  C(N)\ 

Z  Als2  s  S-Pl  S“P2  **•  s“Pn/ 

The  program  first  computes  the  partial  fraction  co¬ 
efficients,  A,  B,  and  C(i)  i  =  1,  2,  . . . ,  N.  The  Z- 
transform  of  each  term  is  then  taken  as  follows: 


(  eg)  )  =  c  u) 

S  -  Pi  Z  -  ePiT 

The  poles  of  G(Z),  exclusive  of  the  one  at  unity,  are 
PiT 

seen  to  be  at  Z  =  e  for  i  =  1,  2,  . ..,  N 

The  program  then  proceeds  to  combine  the  individual 
Z-transforms  over  a  common  denominator,  cancel  the  external 
factor  (1-Z“^),  and  multiply  through  by  the  plant  gain  con¬ 
stant,  VK. 

The  final  results  are  displayed  as  the  coefficients 

of  the  numerator  and  denominator  polynomials  shown  below: 

G(Z)  =  P(PZN  +  P(2)ZN~1  +  P(3)ZN~2  +  +  P(N)Z  +  P(N+1) 

Q(1)ZN+1  +  Q(2)ZN  +  Q( 3)ZN“X  +  ...  +  Q(N+1)Z  +  Q(N+2) 

or,  equivalently: 
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■ft-r  * 


G(Z)  =  PU)^"1  +  P(2)Z~2  +  P(N)Z“N  f  P(N-H)z°(^1} 

Q(1)+Q(2)Z_1 +Q ( 3 ) Z” 2  +  ...  +Q(N+l)Z“N+Q(N+2)Z”tN+1^ 


Discussion  of  the  Data  Provided  by  the  User 

The  general  nature  of  the  input  data  for  each  mode 
is  discussed  below.  The  sequence  in  which  the  data  are 
presented  is  vital ,  and  is  correctly  shown  in  the  discussion. 
Mode  Is 

1)  MODE  =  the  integer  number  1. 

2)  N  =  an  integer  value  between  1  and.  30  representing 
the  degree  of  the  polynomial  whose  roots  are  to  be  found. 

3)  A(i)  =  a  linear  array  of  decimal  values  repre- 
senting  the  coefficients  of  the  polynomial  whose  roots  are 
to  be  found.  These  values  must  appear  on  the  data  card(s) 
from  left  to  right  in  the  same  order  as  they  are  in  the 
polynomial  when  written  in  descending  powers  of  the  vari- 
able.  Example  for  N  =  4j 

A ( 5 ) x4  +  A(4)x3  +  A(3)x2  +  A(2)x  +  A(l)  =  F(x) 

Mode  2s 

1)  MODE  =  the  integer  number  2. 

2)  N  =  an  integer  value  between  1  and  10  representing 
the  degree  of  the  denominator  polynomial  of  the  plant 
transfer  function,,  H(s). 

3)  A(i)  =  a  linear  array  of  decimal  values  repre¬ 
senting  the  coefficients  of  the  denominator  of  H(s)  in 
descending  powers  of  A.  Example  for  N  =  5s 

VK 


H(s)  = 


A(6)s5  +  A(5)s4  +  A(4)s3  +  A( 3 ) s2  +  A(2)s  +  A(l) 
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For  a  type  1  servo,  A(l)  =  0.0  and  N  is  allowed  to 


be  as  high  as  11.  The  pole  at  the  origin  will  be  removed 
by  the  program  and  N  will  be  reduced  by  1  when  a  zero 
value  is  read  for  A(l).  It  is  optional  with  the  user 
whether  or  not  he  removes  the  pole  at  the  origin  himself, 
or  allows  the  program  to  do  it. 

4)  PERIOD  and  VK  =  both  are  decimal  values.  PERIOD 
is  the  reciprocal  of  the  sampling  rate  for  which  the  Z- 
transfer  is  to  be  taken.  VK  is  the  numerator  gain  con¬ 
stant  of  the  plant  transfer  function,  H(s). 

Mode  3j 

1)  MODE  =  the  integer  number  3. 

2)  N  =  the  integer  number  of  poles  of  the  plant 
transfer  function,  exclusive  of  the  pole  at  the  origin. 

3)  P(i)  =  a  linear  array  of  the  decimal  values  of 

the  poles  of  H(s).  Does  not  include  the  pole  at  the  origin. 
All  poles  must  be  simple  and  must  lie  on  the  real  axis  of 
the  s-plane. 

4)  PERIOD  and  VK  =  same  as  in  Mode  2  discussed  above. 

Detailed  Instructions  for  Preparing  Input  Data 

All  input  data  are  provided  to  the  program  through 
the  use  of  IBM  cards.  The  "normal"  number  of  data  cards 
required  by  each  mode  of  operation  is  tabulated  below: 


Mode 

1 

2 

3 


No.  of  data  cards  required 

3 

4 
4 
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"Normally"  it  is  expected  that  the  number  of  decimal 
values  in  the  linear  array  will  be  eight  (8)  or  less,  so 
that  the  entire  array  can  be  contained  on  a  single  card0 
For  arrays  of  more  than  eight  values,  additional  cards 
are  used  as  necessary. 

For  convenience,  the  following  table  summarizes  the 
requirements  for  input  data: 


TABLE  A 


Card  No. 

Datum 

Format 

Used  in 

Mode  1 

Mode  2 

Mode  3 

1 

MODE 

11 

X 

X 

X 

2 

N 

12 

X 

X 

X 

3* 

( Array ) 

8F10.0 

X 

X 

X 

4 

PERIOD,  VK 

F10.0,  F20.0 

X 

X 

*  Use  additional  cards  3  as  necessary  for  arrays  contain¬ 
ing  more  than  eight  (8)  values. 

The  "Format"  information  needs  to  be  explained.  The 
format  letter  I  means  that  the  value  is  an  integer.  In¬ 
tegers  must  be  shown  without  any  decimal  point  and  must 
be  "right  justified"  in  their  allotted  space.  Right 
justified  means  that  the  right-most  digit  of  the  integer 
value  must  appear  in  the  right-most  column  of  the  space 
allotted  for  that  value  on  the  card.  This  applies  only 
to  integer  values. 

The  format  letter  F  indicates  a  decimal  value  and  a 
decimal  point  must  be  used.  If  the  value  is  negative, 
the  minus  sign  must  appear.  Unsigned  values  are  taken 
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to  be  positive.  The  decimal  value,,  including  the  deci¬ 
mal  point  always  and  minus  sign  if  necessary,  may  appear 
anywhere  within  the  allotted  space  (or  "field1"). 

The  number  following  the  letter  I  or  F  is  the  number 
of  columns  on  the  IBM  card  forming  the  allotted  space  or 
field  for  that  value.  There  are  eighty  (80)  columns  on 
a  card  and  they  are  assigned  from  left  to  right  in  order. 

When  a  number  precedes  the  format  letter,  it  tells 
how  many  times  that  particular  field  is  repeated.  Thus 
card  3  contains  eight  identical  fields  of  ten  columns 
each.  The  next  card  may  be  read  by  the  program  as  a  con¬ 
tinuation  of  card  3  if  more  than  eight  values  are  in  the 
array. 

Card  1  uses  only  the  first  column.  This  column  must 
contain  either  the  digit  1,  2,  or  3. 

Card  2  uses  only  the  first  and  second  columns.  If 
the  value  of  N  is  nine  or  less,  then  skip  column  one  and 
place  the  digit  in  column  two. 

Card  3  uses  as  many  identical  fields  of  ten  columns 
each  as  are  necessary.  Note  that  these  fields  begin  in 
columns  1,  11,  21,  31,  etc. 

Card  4  uses  the  first  ten  columns  to  contain  the  deci¬ 
mal  value  for  the  PERIOD.  The  next  twenty  columns  (i.e. , 

11  through  30)  are  reserved  for  the  decimal  value  of  VK. 
Decimal  values  may  lie  anywhere  within  their  allotted 
field  as  long  as  the  decimal  point  is  shown  in  the  field 
and  the  minus  sign  is  in  the  field  for  negative  values. 
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Let  us  now  show  the  preparation  of  data  cards  for 
some  representative  problems. 

EXAMPLE  PROBLEM  -  Mode  1 

Find  all  the  roots  of 

F(x)  =  x7+  16x6+  102x5+  371x4+  838x3+  1200x2+  1012x+  420 
Since  all  we  need  is  a  rootfinding  operation,  we  choose 
Mode  1. 

The  degree,  N,  of  the  polynomial  is  7.  The  polyno- 
mial  is  written  in  descending  powers  of  x  so  the  three 
data  cards  would  be 


Column 

Column 

Column 

Column 

Column 

Column 

Column 

Column 

1 

11 

21 

31 

41 

51 

61 

71 

1 

'1.0 

i 

16.0 

102.0 

371.0 

838.0 

1200.0 

1012.0 

420.0  ( 

Card 
1  3 

J  07 

I 

,  Card 
)  2 

H 

, 

[Card 

1 


The  answers  are: 

Roots  1  and  2  =  -1.0  +  j  1.0 

3  and  4  =  -1.5  +  j  1.658312 

5  and  6  =  -2.0  +  j  1.414214 

7  =  -7.0 


EXAMPLE  PROBLEM  -  Mode  2 


Find  the  Z-transform  of  the  following  configuration 

zero-order 

hold  plant 


O - 

20 

where  the 


■o 


samples/sec  20  samples/sec 

plant  transfer  function  is 
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H(  s )  =  — = -  1,6064235  x_10^ - 

s  +  134s  +  5064s J  +  73274s^  +  356983s 

Since  H(s)  is  not  known  in  factored  form  we  choose 

Mode  2. 

The  degree,  N,  of  the  denominator  of  H(s)  is  5*. 

The  denominator  of  H(s)  is  written  in  descending 
powers.  Note  that  A(l)  is  missing! 

The  numerator  gain  constant,  VK,  is  1.6064235  x  109. 
The  sampling  rate  is  twenty  samples  per  second  so 
the  sampling  PERIOD  is  .05  seconds. 

The  data  cards  are  prepared  as  shown  below? 


Column 

Column 

Column 

Column 

Column  Column 

Column 

Column 

1 

11 

21 

31 

41  51 

61 

71 

* 

1 .05 

1606423500.0 

Card 

4 

>1.0 

i 

134.0 

5064.0 

73274. 

0  356983.0  0.0 

Card 

3 

1  05 

1 

Card 

2 

1  2 
i 

Card 

1 

f 


*  User  may  remove  pole  at  origin,  if  desired,  and  enter  a 

4th  degree  denominator  where  A(l)  =  356983. 

The  answers  are? 

z-transform  =  o(z)  =  ^ 

Q(Z)  Q(1)+Q(2)Z"1+Q(3)Z“^+Q(4)Z“J+ 

P(5)Z~5 

Q(5)Z”4  +  Q(6)Z"5 
where 

P(l)  =  1.54749758  Q(l)  =  1.0 

P( 2 )  =  16.74877866  Q(2)  =  -2.33676594 
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P(3)  =  16.29465286  Q(3)  =  1.92220653 

P( 4 )  =  2.04809669  Q(4)  =  -  .67242310 

P(5)  =  .01866472  Q(5)  =  .08821342 

Q( 6 )  =  -  .00123091 

The  zeros  of  G(Z)  are  at  Z  =  -  .00988230 

-  .13606251 

-  .91925736 
-9.75793506 

The  poles  of  G(Z)  are  at  Z  =  0.57694981 

0.42741498 

0.31663673 

0.01576442 

EXAMPLE  PR03LSM  -  Mode  3 

Find  the  Z-transform  of  the  following  configurations 


zero -order 

hold  plant 


where  the  plant  transfer  function  is 

_  _ 1.6064235  x  109  _ 

'  '  s(s  +  11)  (s  +  17)  (s  +  23)  (s  +  83) 

Since  H(s)  is  known  in  factored  form  we  choose  Mode  3. 

The  number,  N,  of  real  poles  (not  counting  the  one  at 

the  origin)  is  4. 

g 

The  numerator  gain  constant,  VK,  is  1.6064235  x  10  . 

The  sample  PERIOD  is  .05  seconds. 

The  data  cards  are  prepared  as  shown  below: 

Column  Column  Column  Column  Column  Column  Column  Column 

1  11  21  31  41  51  61  71 


i - 

1  .05 

I 


1606423500.0 
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Card 

4 


Column  Column  Column  Column  Column  Column  Column  Column 


1  11  21  31  41  51  61  71 

I - - — - - 

1-11.0  -17.0  -23.0  -83.0  Card 

! _ 3 

1 04  Card 

'  2 

I  1  1  . .  . 

, 3  Card 

i -  1 


The  answers  are  the  same  as  those  obtained  in  the  example 
problem  for  Mode  2. 
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FLOW  DIAGRAM  FOR  PROGRAM  STAPP1 
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APPENDIX  III 


flat j^on  of  stability  of  a  third-order 
System  as  sampling  rate  is  decreased 


or 

P(z) 

Q(») 


cr 

s  +  y+jj-  +”  5 


2-1 

+  70  J 


Ala  4.  Bz 

Fo*  — 


a- 1 


Ca  , 

-  jrr  + 

a-e  a-e 


where  K-  g'O,  000 

As  SI.  /4j?X/0~S~ 
B-  ~3.  yooox  /o~r 
2s_rT  c_-  s.s-ss-t/o-* 

J  2)  =  -0.4O3S*l0~^ 


Mult  iplyi-MCj  out : 

A/  .  3  r-jrr  -?or  _  i  r  -asT  -7or  -nri 

Q(z)_  2  -  [e  +e  f  ljz  f  [e  +e  +-S  J2--<2 


-ff  r 


“7oTi 


P(z)  =  K[AT- Q(iaslecl  i)- c(2+  eOTJ  -  P(^+earTJ]  b-2- 

+  K [- AT (eaf  l eoT) t  &(£ at  +  e7**  i,rr)  +  C ( I  +  ;??dT)  +•  D(i+ Ze**)]  Z 

+  K[ATe9rr-  BeW-  Ce7‘T  -  Pi^T] 


T*(t)  + 


T 


o- 


t-£n 

£0, 00 0 

s 

5’(5>;?.rJ('s>7oy) 

c(t) 
— *»-0 


Characteristic  equation:  A(2)-  P(*)  +  Q(2)  =  £3  *  + 
where  <£?3  =  1 

^  =  A'Ar-  K"  <3  -  KB  ~OT-  JKC-KC  Q°- 3  KD  -  KD~&T- KB  -  e^-^°T-  1 
A 1  =  -KATe*-  KATe°T+  KBZSX+  KB£7ot+  Kb£f’rTi-  KC  +  JKc£~toT 

,  u ,0  .  o  L/n  0“  !s"r  ^*Jrr  -  70 r  -  ?4~r 

+  KD  -r  Or  KD  S  +-  Q  +€  -t-S 

£0  -  KAT e9r-  KB  <s9ST-  KC  e7°-  KD^r-  T 
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Note  that  G(z)  has  been  expressed  in  positive  powers  of 
z  since  we  anticipate  making  a  change  of  variables.  Later 
on„  z~transforms  are  expressed  using  inverse  powers  of  z 
in  order  to  correspond  to  the  concept  of  difference  equa= 
tions. 

_  i.  -h  u> 

Using  the  Bilinear  Transformation?  z  —  i~Z~Zo 

aw  -  &»(■££? +  ^ 

—  $ j  ■*- to  J «/-  [u  40  +  °°  J  f  /  ~ 40  ~  ^  +  10  ] p ~ ~~ 

=  f  to2[345-4t-4+34*]+<o[3a s+ 340]  + 

Routh-Hurwitz  criterion  for  third  order  equations 

all  roots  in  left  half  plane  if  \^02  —  bc  b3  >  O 
to^re  -  &?>  4-  +*  &1  +•  ao 

*b,  =  3^<3  4*  0.%  &  i  “*  3^0 

ba  =  3^3  -tfa  +  3£* 

b3  =  Cl’l  *  H-  j  —  ^  o 
X'K  ter>vt5  <3-f  tlve  5yjte>^  |5Ara-weter53 

fe..  €WV^r;  +i7°T(-KAT)  ♦  k“T(-KAT)  +  (MT-KC-KO) 

bt  =  &r‘(-3KfiT  +  2K&  +  2.)  +  e  (V^T  -  J  KO -3  ) 

t  <£arT(l<AT- 5KS-  3)  4-  (KAT -3KC  -  3KD  -t  a) 

s  ^*'(3  KAT  -4K&-  A)  +  €  (kAT-4*<c) 

4-  e2rT(  KAT  -  4  K D)  4-  (-KAT  tKC+KD  +  l) 
b3  =  e9^T(KAT)  ^'gT(-KAT  UK8t2KC4-2) 

4-  eaST(-KAT  4  3K6  +  2KO  4-  2.)  +  (-  KAT  +-3KC  4  3KD  4  1  ) 

A  computer  program  to  solve  the  above  equations  for 
one  hundred  values  of  T  and  the  results  obtained  are  shown 
on  the  following  pages. 
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. .JOB  STAPP2 

PROGRAM  STAPP2 

C  DETERMINATION  OF  ROUTH-HURWITZ  CRITERION  FOR  THIRD 

ORDER  SYSTEM 

C  AS  SAMPLING  PERIOD  IS  VARIED 

VK=80000. 

A=  . 000571428 
B=-. 000031020 
C-  .000035555 
D=-. 000004535 
T=0 . 0 
N~0 

PRINT  199 

199  FORMAT( 32H  PERIOD  ROUTH  CRITERION) 

100  T~T+. 001 

EXP25s2.71828**(~25.0*T) 

EXP70— 2 , 71828** (-70. 0*T) 

EXP95=2.71828**(-95.0*T) 

AO=EXP95*(VK*A*T-VK*B-l. )-EXP70*VK*C-EXP25*VK*D 
0Al=EXP95*( VK*B+1. )+EXP70*(~VK*A*T+VK*B+2. *VK*C+1. ) 
1  EXP25*(~VK*A*T+VK*B+2.0*VK*D+1. )+VK*C+VK+D 

0A2=EXP70*( -VK*B-VK*C“1.  )+EXP25*(-VK*B~VK*D~l.  )+VK* 
(A*T“2. *C~2. *D) 

1  -1.0-VK*B 

A3— 1. 0 

B0=A3+A2+Al=A0 

Bl=  3 . 0  *A3+A2-Al- 3 . 0  * AO 

B2=3„0*A3‘~A2~Al+3.0*A0 

B3-A3-A2+Al“A0 

ROUTH*  Bl*B2~B0*B3 

PRINT  200, T, ROUTH 

200  FORMAT( F10. 4,  F20.5) 

N=N+1 

IF(N~100 )100,  300, 300 
300  STOP 
END 
END 


(FORTRAN  program  for  the  CDC  1604) 


PERIOD 

ROUTH  CRITERION 

PERIOD 

ROUTH  CRITI 

.0010 

.00060 

.0510 

“2. 33345 

.0020 

.00418 

.0520 

=2.44305 

.0030 

.01226 

.0530 

-2.55283 

.0040 

.02524 

.0540 

“2.66273 

.0050 

.04279 

.0550 

“2.77266 

.0060 

.06411 

.0560 

“2.88256 

.0070 

.08814 

.0570 

“2.99235 

.  0080 

.11372 

.0580 

“3.10199 

.0090 

.13966 

.0590 

=3.21141 

.0100 

.16483 

.0600 

“3.32057 

.0110 

.18818 

.0610 

“3.42940 

.0120 

.20879 

.0620 

“3. 53787 

.0130 

.22587 

.0630 

“3.64593 

.0140 

.23874 

.0640 

=3.75354 

.0150 

.24688 

.0650 

“3.86066 

.0160 

.24986 

.0660 

“3.96725 

.0170 

.24736 

.0670 

“4.07329 

.0180 

.23919 

.0680 

“4.17874 

.0190 

.22522 

.0690 

“4.28358 

.0200 

.20541 

.0700 

“4.38776 

.0210 

.17977 

.0710 

“4.49128 

.0220 

.14839 

.0720 

=4.59410 

.0230 

.11137 

.0730 

“4.69620 

.0240 

.06889 

.0740 

“4.79756 

— ►  .0250 

STABILITY. 02112 

.0750 

“4.89816 

.0260 

LIMIT  ~. 03171 

.0760 

=4.99798 

.0270 

-.08939 

.0770 

-5.09700 

.0280 

-.15167 

.0780 

“5.19522 

.0290 

-.21831 

.0790 

“5.29261 

.0300 

=. 28905 

.0800 

“5.  38915 

.0310 

-.36365 

.0810 

“5.48485 

.0320 

-.44185 

.0820 

“5.57967 

.0330 

=.52341 

.0830 

“5.67363 

.0340 

-.60809 

.0840 

-5.76669 

.0350 

-.69565 

.0850 

“5.85886 

.0360 

“.78588 

.0860 

“5.95012 

.0370 

-.87855 

.0870 

=6.04048 

.0380 

-.97345 

.0880 

“6.12991 

.0390 

“1.07039 

.0890 

“6.21842 

.0400 

“1.16919 

.0900 

=6.30599 

.0410 

-1.26965 

.0910 

“6.39264 

.0420 

“1.37162 

.0920 

“6.47834 

.0430 

“1.47493 

.0930 

“6.56309 

.0440 

“1.57944 

.0940 

“6.64690 

.0450 

-1.68499 

.0950 

“6.72977 

.0460 

=1.79146 

.0960 

“6.81168 

.0470 

“1.89873 

.0970 

“6.89263 

.0480 

“2.00666 

.0980 

“6.97264 

.0490 

=2.11516 

.0990 

“7.05169 

.0500 

-2. 22412 

.1000 

“7.12979 

Determination  of  Routh-Hurwitz  criterion  for  a  third 
order  system  as  sampling  period*  T*  goes  from  .001  to  0.1 
seconds. 
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